Load all important libraries
library(dMod)
library(stringr) # Um bequem mit strings zu arbeiten
library(tidyverse) # Viele Funktionen, u.a. für data.frames und ggplot2 für schöne plots
library(magrittr) # der Pipe-operator %>%: z.B: x = a; y=f(x); z=g(y); wird zu z= a %>% f %>% g
library(conveniencefunctions)
library(libSBML)
load("workspace.rda")
This is basically copied from a libSBML-example
filename <- "/home/denial/Promotion/Projects/methacetin_fitting/model/met13_pkpd_7.xml" # Attention: tilde expansion doesn't work
d <- readSBML(filename)
m <- SBMLDocument_getModel(d)
# errors = SBMLDocument_getNumErrors(d);
# SBMLDocument_printErrors(d);
level = SBase_getLevel (d);
version = SBase_getVersion(d);
cat("File: ",filename," (Level ",level,", version ",version,")\n");
cat(" model id: ", ifelse(Model_isSetId(m), Model_getId(m) ,"(empty)"),"\n");
cat( "functionDefinitions: ", Model_getNumFunctionDefinitions(m) ,"\n" );
cat( " unitDefinitions: ", Model_getNumUnitDefinitions (m) ,"\n" );
cat( " compartmentTypes: ", Model_getNumCompartmentTypes (m) ,"\n" );
cat( " specieTypes: ", Model_getNumSpeciesTypes (m) ,"\n" );
cat( " compartments: ", Model_getNumCompartments (m) ,"\n" );
cat( " species: ", Model_getNumSpecies (m) ,"\n" );
cat( " parameters: ", Model_getNumParameters (m) ,"\n" );
cat( " initialAssignments: ", Model_getNumInitialAssignments (m) ,"\n" );
cat( " rules: ", Model_getNumRules (m) ,"\n" );
cat( " constraints: ", Model_getNumConstraints (m) ,"\n" );
cat( " reactions: ", Model_getNumReactions (m) ,"\n" );
cat( " events: ", Model_getNumEvents (m) ,"\n" );
cat( "\n" );
# Mit dem Pipe-Operator kann man Funktionen verketten
# Standardmäßig wird das vorherige Ergebnis als erstes Argument von der nächsten Funktion eingesetzt. wenn man das nicht will, kann man es als . woanders hinsetzen
f <- function(x) x^2;
g <- function(x,y) x-y;
2 %>% f # f(2)
2 %>% f %>% g(3) # g(f(2),3) = 4-3 = 1
2 %>% f %>% g(3,.) # g(3,f(2)) = 3-4 = -1
# Man kann auch Funktionen definieren, die mit . losgehen, was dann das Argument für die Funktion ist.
h <- . %>% sqrt %>% add(5)
h(4)
# Das ist besonders nützlich in lapply, sapply und so weiter, wo man über eine liste(oder einen vektor) immer wieder die gleiche funktion laufen lässt (syntactic sugar für for-loops)
sapply(1:4, function(i) i^2 +5)
sapply(1:4 , . %>% raise_to_power(2) %>% add(5) ) # dassselbe
# get rules
nrules <- Model_getNumRules(m)
lrules <- Model_getListOfRules(m)
rules <- structure(sapply(0:(nrules-1), . %>% ListOfRules_get(lrules,.) %>% Rule_getFormula),
names = sapply(0:(nrules-1), . %>% ListOfRules_get(lrules,.) %>% Rule_getId))
rulenames <- names(rules)
# "Cure" rules: Since I do a parameter trafo for the units, I don't want to have any unit conversions via some rules
# A "bad" rule would be eg "QC = CO*3600/100", since I take care of the Units later. Therefore, the rule should be only "QC = CO"
rules <- rules %>% str_replace_all(c("1000" = "1", "3600" = "1", "\\b60\\b" = 1)) %>% set_names(rulenames)
# Apply the rules onto themselves to insert parameter transformations
# Final goal is to have a named vector where
# names are the "inner" parameters that are used within the model
# values are functions of "outer" parameters that are fed into the model
# apply rules 1st time
rules <- paste0("(", rules, ")") %>% set_names(paste0("\\b", rulenames, "\\b"))
rules <- str_replace_all(rules, rules) %>% set_names(rulenames)
# apply rules 2nd time
rules <- paste0("(", rules, ")") %>% set_names(paste0("\\b", rulenames, "\\b"))
rules <- str_replace_all(rules, rules) %>% set_names(rulenames)
# check if any of the rules are functions of other rules
indices <- rules %>% sapply(. %>% str_detect(paste0("\\b",rulenames, "\\b")) %>% any)
rules %>% extract(indices) %>% sapply(. %>% str_detect(paste0("\\b",rulenames, "\\b")) %>% extract(rulenames,.)) # check works, none of the
# print(rules)
# getSymbols(rules) # These are the "outer" parameters
# get reactions
nreactions <- Model_getNumReactions(m)
lreactions <- Model_getListOfReactions(m)
reactions <- lapply(0:(nreactions-1), function(i) {
reaction <- ListOfReactions_get(lreactions,i)
nreactants <- reaction %>% Reaction_getNumReactants()
if (nreactants > 0) {
lreactants <- reaction %>% Reaction_getListOfReactants()
myreactant_stoichiometries <- lapply(0:(nreactants-1), . %>% ListOfSpeciesReferences_get(lreactants,.) %>% SpeciesReference_getStoichiometry)
myreactant_IDs <- lapply(0:(nreactants-1), . %>% ListOfSpeciesReferences_get(lreactants,.) %>% SimpleSpeciesReference_getSpecies)
from = paste0(paste0("(", myreactant_stoichiometries, "*", myreactant_IDs, ")"), collapse = "+")
} else {
from = ""
}
nproducts <- reaction %>% Reaction_getNumProducts()
if (nproducts > 0) {
lproducts <- reaction %>% Reaction_getListOfProducts()
myproduct_stoichiometries <- lapply(0:(nproducts-1), . %>% ListOfSpeciesReferences_get(lproducts,.) %>% SpeciesReference_getStoichiometry)
myproduct_IDs <- lapply(0:(nproducts-1), . %>% ListOfSpeciesReferences_get(lproducts,.) %>% SimpleSpeciesReference_getSpecies)
to = paste0(paste0("(", myproduct_stoichiometries, "*", myproduct_IDs, ")"), collapse = "+")
} else {
to = ""
}
# Apply rules to rate expressions
myrules <- rules
mynames <- names(myrules)
absorption_indices <- str_detect(myrules, "Absorption")
absorption_rules <- myrules[absorption_indices] %>% structure(names = mynames[absorption_indices])
myrules <- structure(myrules[!str_detect(myrules, "Absorption")], names = paste0("\\b", mynames[!str_detect(myrules, "Absorption")], "\\b"))
rate <- reaction %>% Reaction_getKineticLaw() %>% KineticLaw_getFormula() %>% str_replace_all(myrules)
description <- reaction %>% Reaction_getId()
# Incorporate the absorption:
# For this I add another reaction which absorbs the oral dose, e.g. D_apap_sul
if(description %>% str_detect("Absorption")) {
from[2] <- names(absorption_rules)[absorption_rules %>% str_detect(paste0(description, "\\)*$"))]
to[2] <- ""
rate[2] <- rate
description[2] <- description
}
# print(i)
return(data.frame(from = from, to = to, rate = rate, description = description, stringsAsFactors = F))
}) %>% do.call(rbind,.) # make one big data.frame out of it
# Build the dMod-object "eqnlist", which stores the reactions and stoichiometries
el <- eqnlist()
for(i in 1:nrow(reactions)) el <- addReaction(el, reactions$from[i], reactions$to[i], reactions$rate[i], reactions$description[i])
# Convert to "eqnvec", which is basically a named vector of the ODEs and the names denote the states
f <- el %>% as.eqnvec()
# This is the ODE when every rule is applied
print(f %>% head)
# get the parameters from the definition
nsbml_parameters <- m %>% Model_getNumParameters()
lsbml_parameters <- m %>% Model_getListOfParameters()
sbml_parameters <- structure(sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getValue),
names = sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getId))
sbml_parameter_descriptions <- structure(sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getName),
names = sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getId)) %>% {.[order(names(.))]}
# unit conversion
sbml_parameter_units <- structure(sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getUnits), names = sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getId))
# The factors to bring each parameter to the units seconds, grams, litres
multiplication_factors <- sapply(0:(nsbml_parameters-1), . %>%
ListOfParameters_get(lsbml_parameters,.) %>%
Parameter_getDerivedUnitDefinition %>%
{myunitdef <- .
nunits <- UnitDefinition_getNumUnits(myunitdef)
lunits <- UnitDefinition_getListOfUnits(myunitdef)
lapply(0:(nunits-1), . %>%
ListOfUnits_get(lunits,.) %>%
{(Unit_getMultiplier(.) * 10^(Unit_getScale(.)))^Unit_getExponent(.)}) %>% Reduce("*",.)
} %>%
{.}
)
sbml_parameters <- sbml_parameters * multiplication_factors
# Initial assignments
niass <- m %>% Model_getNumInitialAssignments()
liass <- m %>% Model_getListOfInitialAssignments()
iass <- structure(sapply(0:(niass-1), . %>% ListOfInitialAssignments_get(liass,.) %>% print %>% str_split("\n") %>% sapply( . %>% {.[str_detect(.,"<ci>")]}) %>% str_replace_all("^.*<ci> ", "") %>% str_replace_all(" </ci>.*$", "")), names = sapply(0:(niass-1), . %>% ListOfInitialAssignments_get(liass,.) %>% InitialAssignment_getId))
nspecies <- Model_getNumSpecies(m)
lspecies <- Model_getListOfSpecies(m)
species <- lapply(0:(nspecies-1), . %>% {ListOfSpecies_get(lspecies,.)} %>% Species_getId()) %>% do.call(c,.)
inits <- structure(sapply(0:(nspecies-1), . %>% ListOfSpecies_get(lspecies,.) %>% Species_getInitialAmount), names = species)
inits[names(iass)] <- sbml_parameters[iass]
# All parameters combined, more than actually needed, because many of them have been replaced due to Rules
all_pars <- c(sbml_parameters, inits) # all possible parameters
# To check if the unit conversion worked
data.frame(par = names(sbml_parameters), value = sbml_parameters, original_unit= sbml_parameter_units, multiplication_factor = multiplication_factors, original_value = sbml_parameters/multiplication_factors, row.names = NULL) %>%
head()
#%>%
# filter(par %in% c("MET2APAP_HLM_CL", "fumic_metc13", "MPPGL", "BW", "FVli")) #%>% summarise(prod(multiplication_factor))
# filter(par %in% c("CO", "QC"))
# rules[str_detect(rules, "MET2APAP")]
# rules[str_detect(rules, "CO")]
the odemodel() command takes a while because sensitivity equations are calculated for derivatives and then the whole system is compiled into c-code.
# myodemodel <- odemodel(f, modelname = "methacetin") # this compiles the ode into a c-file. you can comment out this and the next line after it has been run once to save time.
# save(myodemodel, file = "methacetin.rda")
load("methacetin.rda")
Make a prediction function from the odemodel. x will be a function x(times, pars)
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
This plot is supposed to be the first plot from chunk 3 of the html-file that you sent me via email.
pars["Ave_metc13"] <- 0
pars["D_metc13"] <- 0.1 #100mg
mytimes = seq(0,8*3600, length.out = 80) #8 hours
# pred <- (g*x)(0:50, mypars, deriv = F, conditions = rownames(attr(data, "cond")))
prediction <- (x)(times = mytimes, pars = pars, deriv = F, conditions = "c1")
# Compute the volumes in litres for each state
volumes <- lapply(0:(nspecies-1), . %>% {ListOfSpecies_get(lspecies,.)} %>% Species_getCompartment()) %>% do.call(c,.) %>% set_names(species)
myrules <- rules %>% set_names(paste0("\\b", rulenames, "\\b"))
volumes <- str_replace_all(volumes, myrules) %>% set_names(species)
vol_fun <- funC0(volumes)
vol <- do.call(vol_fun, as.list(pars)) %>% t %>% {data.frame(volume=., name = rownames(.), stringsAsFactors = F)}
# Plot
plot_fun <- function(pred) pred %>%
wide2long() %>%
full_join(vol,by= "name") %>%
mutate(value = value * 1000, time = time/3600, concentration = value/volume) %>% # scale g to mg, s to h
separate(col = name, into = c("compartment", "substance"), sep = "_", extra = "merge") %>%
filter(compartment %in% c("Ali", "Ave", "Alu")) %>% # plot only liver, venuous and lung (as in the html file)
ggplot(aes(x= time, y = concentration)) +
geom_line(aes(color = compartment)) +
facet_wrap(~substance, scales = "free") +
theme_dMod() + scale_color_dMod()
plot_fun(prediction)
prediction_long <- (x)(times = mytimes*20,pars = pars, deriv = F, conditions = "c1")
plot_fun(prediction_long)
# Second plot from chunk 3
pars <- all_pars[getParameters(x)]
pars["Ave_metc13"] <- 0.1
pars["D_metc13"] <- 0 #100mg
# pars["Kplu_metc13"] <- 1
# pars["FVlu"] <- 7.6*10^(-6) * 10
prediction_second_plot <- (x)(times = mytimes,pars = pars, deriv = F, conditions = "c1")
plot_fun(prediction_second_plot)
Reduce the model complexity by inserting the fixed parameter values
free_parameters <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km" # Km value
)
fixed_parameters <- pars[!(names(pars)%in%c(free_parameters,names(f)[1]))] %>% names
myodemodel <- odemodel(f, modelname = "methacetin_reduced", fixed = fixed_parameters)
x <- Xs(myodemodel)
observables <- c(apap = "Ave_apap/(BW*FVve)",
apap_glu = "Ave_apap_glu/(BW*FVve)",
apap_sul = "Ave_apap_sul/(BW*FVve)",
apap_cys = "Ave_apap_cys/(BW*FVve)"
)
g <- Y(observables, x, parameters = free_parameters)
(g*x)(mytimes, pars) %>% set_names("condition1") %>% plotPrediction(name %in% names(observables))
Load the data and transform it
myfiles <- list.files("~/Promotion/Projects/methacetin_fitting/data/", full.names = T)
raw_data <- myfiles %>% lapply(. %>% read.table(header = T, sep = "\t", stringsAsFactors = F))
data <-
raw_data %>%
lapply(. %>%
select(-contains("_mol")) %>%
gather("name_std", "std", ends_with("_sd")) %>%
mutate(name_std = str_replace(name_std, "_sd","")) %>%
gather("name_sigma", "sigma", ends_with("_se")) %>%
mutate(name_sigma = str_replace(name_sigma, "_se","")) %>%
{gather(.,"name", "value", one_of(.$name_std))} %>%
filter(name == name_std, name == name_sigma) %>%
# select(name,time,value,sigma) %>%
{.}) %>%
do.call(dMod::combine,.) %>%
mutate(D_apap = "D_apap", Ave_apap = "Ave_apap" ) %>%
{.$D_apap[.$study=="Chan1997"] <- 1400 / 1000
.$Ave_apap[.$study=="Chan1997"] <- 0
.$D_apap[.$study=="Chiew2010"] <- 5600 / 1000
.$Ave_apap[.$study=="Chiew2010"] <- 0
.$D_apap[.$study=="Critchley2005"] <- 1400 /1000
.$Ave_apap[.$study=="Critchley2005"] <- 0
.$D_apap[.$study=="Rawlins1977"] <- .$dose[.$study=="Rawlins1977"] * (.$route[.$study=="Rawlins1977"]%in%"oral") / 1000
.$Ave_apap[.$study=="Rawlins1977"] <- .$dose[.$study=="Rawlins1977"] * (.$route[.$study=="Rawlins1977"]%in%"iv") / 1000
.
} %>%
mutate(time = time * 3600, value = value/1000, sigma = sigma/1000) %>%
select(-group, -health_status, - name_std, - name_sigma, -std, -ethnicity, -route, -dose, -substance) %>%
# filter(!is.na(sigma)) %>%
# as.datalist() %>%
{.}
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
plot(mydatalist)
Parameter transformations to define the conditions
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters] <- paste0("exp(log", free_parameters, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p <- NULL
for(i in 1:length(p_list)) { p <<- p + p_list[[i]]}
pouter <- log(pars[free_parameters]) %>% set_names(paste0("log",names(.)))
mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
plotCombined(mypred, mydatalist, name%in% names(observables))
obj <- normL2(mydatalist, (g*x*p))
# myfit <- mstrust(objfun = obj, center = pouter, studyname = "methacetin", cores = 3)
# save(myfit, file = "fit.rda")
load("fit.rda")
fitted_pars <- myfit %>% as.parframe() %>% as.parvec()
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
mypred <- (g*x*p)(seq(0, 48*3600, length.out = 100), fitted_pars)
liver <- names(f)[str_detect(names(f), "li")]
medication <- c(names(f)[str_detect(names(f), "D_apap$")], "Ave_apap")
myfit %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
myfit %>% as.parframe() %>% plotValues()
plotCombined(mypred, mydatalist, name %in% c(names(observables), liver, medication))
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 8 2407.821 TRUE 60 -5.125968 -2.129845
## 2 35 2407.821 TRUE 71 -5.125973 -2.129851
## 3 21 2407.821 TRUE 39 -5.125965 -2.129842
## 4 29 2407.821 TRUE 67 -5.125962 -2.129839
## 5 33 2407.821 TRUE 76 -5.125966 -2.129843
## 6 22 2407.821 TRUE 75 -5.125969 -2.129847
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 -6.118641 -10.49793 -6.301154
## 2 -6.118641 -10.49793 -6.301161
## 3 -6.118641 -10.49794 -6.301223
## 4 -6.118641 -10.49794 -6.301221
## 5 -6.118641 -10.49794 -6.301214
## 6 -6.118641 -10.49794 -6.301198
## logAPAPCYS_HLM_CL: 2.75936204923354e-05
## logAPAPCYS_Km: 0.00183418625326802
## logAPAPGLU_HLM_CL: 0.00594046463001759
## logAPAPGLU_Km: 0.118855661904664
## logAPAPSUL_HLM_CL: 0.00220144570190944
Here I freed some additional parameters, but since I freed “Ka_apap_sul” instead of “Ka_apap”, the results don’t make much additional sense.
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters1 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap_sul", "F_apap_sul",
"Kpre_apap", "Kpki_apap", "Kpli_apap",
"Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
"Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
"Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu",
"Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
"Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters1 <- pars[!(names(pars)%in%c(free_parameters1,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters1] <- paste0("exp(log", free_parameters1, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p1 <- NULL
for(i in 1:length(p_list)) { p1 <<- p1 + p_list[[i]]}
pouter <- log(pars[free_parameters1]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj1 <- normL2(mydatalist, (g*x*p1))
# job1 <- runbg({myfit <- mstrust(objfun = obj1, center = pouter, studyname = "methacetin", cores = 20, fits = 100); myfit}, machine = "knecht3", filename = "job1")
# myfit1 <- job1$get()$knecht3
# save(myfit1, file = "myfit1.rda")
# job1$purge()
load("myfit1.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit1 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit1 %>% as.parframe())+scale_y_log10()
mypred1 <- (g*x*p1)(mytimes*4, myfit1 %>% as.parframe() %>% as.parvec)
plotCombined(mypred1, mydatalist, name %in% names(observables))
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 71 2749.704 FALSE 100 -3.121425 0.9276677
## 2 17 2871.109 FALSE 100 -4.825791 -0.8966849
## 3 4 2879.499 TRUE 82 -4.836514 -2.0500243
## 4 83 3006.360 TRUE 38 -4.810375 -1.6071039
## 5 11 3025.134 TRUE 47 -5.068367 -2.1571283
## 6 60 3063.915 FALSE 100 -4.778175 -1.9189665
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1 -5.966762 -8.877054 -8.056205 -9.043602
## 2 -7.312292 -9.776533 -6.317704 -7.411144
## 3 -6.424722 -9.745092 -7.247758 -7.449345
## 4 -6.432805 -10.767499 -6.451981 -6.734894
## 5 -6.368207 -9.236633 -7.298571 -8.975554
## 6 -6.479809 -10.151009 -5.692777 -5.931667
## logF_apap_sul logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys
## 1 0.53130523 -0.2895295 -0.9819923 0.32392147 3.4721895
## 2 0.07956807 -0.3510074 0.2160174 0.99467641 0.7860616
## 3 -0.52122381 -0.3932631 -1.0615523 -0.09235355 0.9514369
## 4 0.16252315 -0.3785620 0.1383524 0.27175050 -0.8202826
## 5 1.21332440 -0.2100422 -1.7021606 -0.01216570 1.9565013
## 6 -2.38273975 -0.4474130 -1.6230284 0.04538435 0.4837584
## logKpki_apap_cys logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu
## 1 -0.9981306 -0.7078988 -1.1713447 -0.4586217
## 2 1.4876644 0.5391938 0.5938409 0.6894229
## 3 1.3055012 0.4413984 0.7986194 0.3810163
## 4 -0.1222928 -0.2784549 -0.6756876 0.1384249
## 5 1.2069110 0.6198618 0.5661661 -1.1336736
## 6 -0.1000206 -0.1492559 0.6852967 -0.2904269
## logKpli_apap_glu logKpre_apap_sul logKpre_co2c13 logKpli_co2c13
## 1 1.2820601 0.8610661 -1.8050599 -1.0958784
## 2 -1.2743622 -2.0279402 -1.6334736 0.2964431
## 3 -0.7746656 -1.5237195 0.3130873 0.3634344
## 4 0.4949366 -0.5271704 -0.8464643 -1.1150170
## 5 -1.4231745 -1.4809556 0.9901271 -0.4769397
## 6 1.7057846 -1.3223539 -1.2621168 -0.9185781
## logKpre_metc13 logKpli_metc13
## 1 -2.4242875 -0.4355059
## 2 -1.3074538 -0.5878608
## 3 0.0672823 1.2834798
## 4 -2.5253310 1.1213766
## 5 -1.0586655 1.3624494
## 6 -1.3302118 1.4817382
## logAPAPCYS_HLM_CL: 0.000139554700132717
## logAPAPCYS_Km: 0.000317128135451512
## logAPAPGLU_HLM_CL: 0.0440942903336689
## logAPAPGLU_Km: 2.52860494207782
## logAPAPSUL_HLM_CL: 0.00256252477737589
## logF_apap_sul: 1.70115125864559
## logKa_apap_sul: 0.000118144550067709
## logKpki_apap: 0.374564129323335
## logKpki_apap_cys: 0.368567801198445
## logKpki_apap_glu: 0.63215432668864
## logKpli_apap: 1.38253873675798
## logKpli_apap_cys: 0.492678343277563
## logKpli_apap_glu: 3.60405672997012
## logKpli_co2c13: 0.334245860035197
## logKpli_metc13: 0.646937302760331
## logKpre_apap: 0.748615739138366
## logKpre_apap_cys: 32.2071824118807
## logKpre_apap_glu: 0.309949884477503
## logKpre_apap_sul: 2.36568138653131
## logKpre_co2c13: 0.164464609213266
## logKpre_metc13: 0.0885411804880305
Here I freed some less parameters than in try 1 but it suffers from the same problem
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters2 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap_sul", "F_apap_sul" #, # total dumm, Ka und F sind so redundant wie es nur geht.
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu",
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters2 <- pars[!(names(pars)%in%c(free_parameters2,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters2] <- paste0("exp(log", free_parameters2, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p2 <- NULL
for(i in 1:length(p_list)) { p2 <<- p2 + p_list[[i]]}
pouter <- log(pars[free_parameters2]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj2 <- normL2(mydatalist, (g*x*p2))
# job2 <- runbg({myfit <- mstrust(objfun = obj2, center = pouter, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht4", filename = "job2")
#
# myfit2 <- job2$get()$knecht4
# save(myfit2, file = "myfit2.rda")
# job2$purge()
load("myfit2.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit2 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit2 %>% as.parframe())+scale_y_log10()
mypred2 <- (g*x*p2)(mytimes*4, myfit2 %>% as.parframe() %>% as.parvec)
plotCombined(mypred2, mydatalist, name %in% names(observables))
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 24 3134.039 TRUE 25 -5.067594 -2.064909
## 2 56 3134.040 TRUE 30 -5.067589 -2.064904
## 3 20 3134.040 TRUE 26 -5.067595 -2.064908
## 4 89 3134.040 TRUE 16 -5.067592 -2.064905
## 5 37 3134.040 TRUE 86 -5.067589 -2.064901
## 6 32 3134.040 TRUE 20 -5.067590 -2.064902
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1 -6.117757 -10.48519 -6.241245 -7.346186
## 2 -6.117757 -10.48520 -6.241296 -7.021465
## 3 -6.117758 -10.48501 -6.240447 -7.216762
## 4 -6.117758 -10.48501 -6.240447 -6.306130
## 5 -6.117757 -10.48501 -6.240481 -7.986534
## 6 -6.117757 -10.48501 -6.240461 -7.899805
## logF_apap_sul
## 1 41.111322
## 2 43.933670
## 3 50.970546
## 4 43.726021
## 5 7.397659
## 6 34.039840
## logAPAPCYS_HLM_CL: 2.79472212506865e-05
## logAPAPCYS_Km: 0.00194743024961977
## logAPAPGLU_HLM_CL: 0.00629755070545838
## logAPAPGLU_Km: 0.126829778557011
## logAPAPSUL_HLM_CL: 0.00220339319347716
## logF_apap_sul: 715187883400265600
## logKa_apap_sul: 0.000645048206917656
Same problem as in tries 1 and 2
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters3 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap_sul", "F_apap_sul" ,
"Kpre_apap", "Kpki_apap", "Kpli_apap",
"Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
"Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
"Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters3 <- pars[!(names(pars)%in%c(free_parameters3,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters3] <- paste0("exp(log", free_parameters3, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p3 <- NULL
for(i in 1:length(p_list)) { p3 <<- p3 + p_list[[i]]}
pouter <- log(pars[free_parameters3]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj3 <- normL2(mydatalist, (g*x*p3))
# job3 <- runbg({myfit <- mstrust(objfun = obj3, center = pouter, studyname = "methacetin", cores = 15, fits = 100); myfit}, machine = "knecht1", filename = "job3")
# myfit3 <- job3$get()$knecht1
# save(myfit3, file = "myfit3.rda")
# job3$purge()
load("myfit3.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit3 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit3 %>% as.parframe())+scale_y_log10()
mypred3 <- (g*x*p3)(mytimes*4, myfit3 %>% as.parframe() %>% as.parvec)
plotCombined(mypred3, mydatalist, name %in% names(observables))
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 59 2562.482 TRUE 84 -5.201098 -2.097660
## 2 60 2689.397 TRUE 48 -5.383774 -2.043510
## 3 55 2856.425 FALSE 100 -5.109883 -1.468816
## 4 48 3019.107 FALSE 100 -6.865873 -5.295360
## 5 1 3039.058 FALSE 100 -4.638935 -1.383292
## 6 89 3066.281 TRUE 76 -6.464966 -4.620926
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1 -6.555927 -10.709675 -6.33769669 -6.042436
## 2 -6.862922 -10.910624 -6.20766738 -8.289160
## 3 -6.656348 -9.872593 -4.75970985 -6.579562
## 4 -6.050377 -9.618928 -6.86802281 -6.542891
## 5 -6.709802 -7.521883 -0.09424782 -7.980098
## 6 -5.857182 -10.682944 -8.86895808 -7.943206
## logF_apap_sul logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys
## 1 2.808457812 -0.3669234 0.7613894 0.1672635 -0.6048195
## 2 -1.512335536 -0.1336978 -0.2052586 0.5449920 -1.5425344
## 3 0.009716211 -0.2655233 -1.6898116 0.6768234 -0.2769801
## 4 -0.661252544 -0.1334166 0.2229806 -0.6106942 -0.4045914
## 5 2.748915536 -0.4468966 -0.6003575 0.3837791 -0.3193860
## 6 1.679453017 -0.1192626 -0.5590351 -0.6042076 -2.9191368
## logKpki_apap_cys logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu
## 1 -0.1241948 -0.04614371 -0.8681152 0.36911867
## 2 -0.2511150 1.52077842 -0.8083817 0.16712529
## 3 0.7221217 -0.09379088 -1.0941684 0.26637337
## 4 2.2752818 0.30050894 -1.5193795 1.64286021
## 5 -2.5272957 1.04456393 0.3086699 0.04221941
## 6 0.4521719 -0.72613317 0.1994445 0.85658306
## logKpli_apap_glu logKpre_apap_sul
## 1 -0.003140001 -1.0894804
## 2 -0.342208285 -2.2243255
## 3 -0.207101279 -0.1223363
## 4 -0.702140533 -2.0416762
## 5 0.408290615 -1.1143709
## 6 -0.820070493 -0.7383285
## logAPAPCYS_HLM_CL: 2.23278620657747e-05
## logAPAPCYS_Km: 0.00176837065068179
## logAPAPGLU_HLM_CL: 0.00551051224903083
## logAPAPGLU_Km: 0.122743302983146
## logAPAPSUL_HLM_CL: 0.0014216640114061
## logF_apap_sul: 16.5843223443174
## logKa_apap_sul: 0.00237576415394654
## logKpki_apap: 2.14124916301601
## logKpki_apap_cys: 0.883207793863882
## logKpki_apap_glu: 1.44645924035911
## logKpli_apap: 1.18206573439193
## logKpli_apap_cys: 0.954904720644904
## logKpli_apap_glu: 0.996864923604179
## logKpre_apap: 0.692862701926673
## logKpre_apap_cys: 0.546172985593896
## logKpre_apap_glu: 0.419741938691713
## logKpre_apap_sul: 0.336391236317653
Same problem as in tries 1,2,3. The only difference to try 2 is that I removed the structural non-identifiability of Ka_apap_sul and F_apap_sul.
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters4 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap_sul"#, "F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters4 <- pars[!(names(pars)%in%c(free_parameters4,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters4] <- paste0("exp(log", free_parameters4, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p4 <- NULL
for(i in 1:length(p_list)) { p4 <<- p4 + p_list[[i]]}
pouter <- log(pars[free_parameters4]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj4 <- normL2(mydatalist, (g*x*p4))
# job4 <- runbg({myfit <- mstrust(objfun = obj4, center = pouter, studyname = "methacetin", cores = 15, fits = 100); myfit}, machine = "knecht1", filename = "job4")
# save(job4, file = "job4.rda")
# myfit4 <- job4$get()$knecht1
# save(myfit4, file = "myfit4.rda")
# job4$purge()
load("myfit4.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit4 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit4 %>% as.parframe())+scale_y_log10()
mypred4 <- (g*x*p4)(mytimes*4, myfit4 %>% as.parframe() %>% as.parvec, deriv = F)
myplot4 <- plotCombined(mypred4, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot4)
# myplot4
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 65 3134.039 TRUE 31 -5.067596 -2.064909
## 2 47 3134.040 TRUE 74 -5.067589 -2.064901
## 3 51 3134.040 TRUE 82 -5.067591 -2.064903
## 4 49 3134.040 TRUE 58 -5.067589 -2.064902
## 5 17 3134.040 TRUE 18 -5.067592 -2.064904
## 6 91 3134.040 TRUE 25 -5.067590 -2.064903
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1 -6.117758 -10.48501 -6.240449 69.0940136
## 2 -6.117757 -10.48502 -6.240488 -1.9728623
## 3 -6.117757 -10.48501 -6.240450 -1.9723416
## 4 -6.117757 -10.48501 -6.240446 -0.7434353
## 5 -6.117757 -10.48501 -6.240455 21.5604671
## 6 -6.117758 -10.48499 -6.240395 39.2232311
## logAPAPCYS_HLM_CL: 2.79524088116887e-05
## logAPAPCYS_Km: 0.00194898079879723
## logAPAPGLU_HLM_CL: 0.0062975431281648
## logAPAPGLU_Km: 0.126829823077826
## logAPAPSUL_HLM_CL: 0.00220339066957284
## logKa_apap_sul: 1.01659705294599e+30
Freeing the parameters in the upper four fits didn’t help much. On the other hand, lots of the fits didn’t converge. From here, I could do several things.
Here I freed Ka_apap which makes sense conceptually (increase absorption rate). But it doesn’t improve the fits much.
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters5 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, "F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters5 <- pars[!(names(pars)%in%c(free_parameters5,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters5] <- paste0("exp(log", free_parameters5, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p5 <- NULL
for(i in 1:length(p_list)) { p5 <<- p5 + p_list[[i]]}
pouter <- log(pars[free_parameters5]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj5 <- normL2(mydatalist, (g*x*p5))
# job5 <- runbg({myfit <- mstrust(objfun = obj5, center = pouter, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht1", filename = "job5")
# save(job5, file = "job5.rda")
# myfit5 <- job5$get()$knecht1
# save(myfit5, file = "myfit5.rda")
# job5$purge()
load("myfit5.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit5 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit5 %>% as.parframe())+scale_y_log10()
mypred5 <- (g*x*p5)(mytimes*5, myfit5 %>% as.parframe() %>% as.parvec, deriv = F)
myplot5 <- plotCombined(mypred5, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot5)
# myplot5
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 33 3092.015 FALSE 100 -5.048048 -2.013030
## 2 12 3092.025 TRUE 51 -5.047014 -2.013058
## 3 99 3092.031 FALSE 100 -5.048425 -2.014686
## 4 60 3092.117 TRUE 96 -5.047379 -2.011610
## 5 88 3092.146 TRUE 49 -5.047300 -2.011395
## 6 26 3092.203 TRUE 40 -5.048943 -2.016184
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -6.142913 -10.48392 -6.226134 -6.853661
## 2 -6.141722 -10.48302 -6.222827 -6.879080
## 3 -6.141694 -10.48360 -6.225265 -6.880146
## 4 -6.143583 -10.48343 -6.224064 -6.838579
## 5 -6.143710 -10.48332 -6.223574 -6.835689
## 6 -6.140809 -10.48370 -6.225923 -6.898840
## logAPAPCYS_HLM_CL: 2.79827234659552e-05
## logAPAPCYS_Km: 0.00197708080993727
## logAPAPGLU_HLM_CL: 0.00642185876478594
## logAPAPGLU_Km: 0.133583353550372
## logAPAPSUL_HLM_CL: 0.0021486564154929
## logKa_apap: 0.00105558386269278
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters6 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap", #"F_apap_sul" ,
"Kpre_apap", "Kpki_apap", "Kpli_apap",
"Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
"Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
"Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters6 <- pars[!(names(pars)%in%c(free_parameters6,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters6] <- paste0("exp(log", free_parameters6, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p6 <- NULL
for(i in 1:length(p_list)) { p6 <<- p6 + p_list[[i]]}
pouter <- log(pars[free_parameters6]) %>% set_names(paste0("log",names(.)))
best_fit <- myfit %>% as.parframe() %>% as.parvec()
pouter[names(best_fit)] <- best_fit
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj6 <- normL2(mydatalist, (g*x*p6))
# job6 <- runbg({myfit <- mstrust(objfun = obj6, center = pouter, studyname = "methacetin", cores = 10, fits = 100); myfit}, machine = "knecht3", filename = "job6")
# save(job6, file = "job6.rda")
# job6$check()
# myfit6 <- job6$get()$knecht3
# save(myfit6, file = "myfit6.rda")
# job6$purge()
load("myfit6.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit6 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit6 %>% as.parframe())+scale_y_log10()
mypred6 <- (g*x*p6)(mytimes*4, myfit6 %>% as.parframe() %>% {.[2,]} %>% as.parvec)
myplot6 <- plotCombined(mypred6, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot6)
# myplot6
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 16 2654.322 TRUE 88 -6.288191 -4.0064053
## 2 60 3025.841 FALSE 100 5.036281 8.1988600
## 3 39 3059.748 FALSE 100 -3.427671 0.2086894
## 4 19 3210.587 TRUE 100 -4.251426 -1.6417543
## 5 47 3233.777 TRUE 48 -5.809100 -2.7455612
## 6 48 3302.297 FALSE 100 -4.770855 -0.8605885
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -6.242811 -9.413691 -6.941120 -5.245524
## 2 -5.281948 -10.686351 -6.408999 -7.197151
## 3 -6.146221 -9.736307 -7.455677 -6.582344
## 4 -4.910387 -10.540902 -7.410143 -7.443229
## 5 -6.315936 -10.030452 -1.701016 -7.739147
## 6 -7.215143 -9.338544 -1.877134 -7.204190
## logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys logKpki_apap_cys
## 1 -0.04690208 0.1358953 -0.3757922 2.19032710 -0.4659840
## 2 -0.27191200 1.4639032 -0.5538510 -0.52221741 -0.3912030
## 3 -0.42939025 0.3304541 0.2229900 0.05838160 1.8274448
## 4 -0.18759530 -1.7257917 -0.7659503 -0.40658360 0.1236754
## 5 -0.08418012 0.4028411 0.2491479 -3.17511619 -2.8534696
## 6 -0.34371554 -1.2734140 1.0683678 0.03096511 -1.5696406
## logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu logKpli_apap_glu
## 1 -0.298471419 -1.2802038 0.54299839 0.09378855
## 2 0.276192230 -1.6506727 -0.45872701 -0.11353331
## 3 0.279366410 -1.3550508 -0.31509066 -0.10951404
## 4 0.009838299 -0.8169902 -0.19239433 -0.32722341
## 5 -0.007563093 -2.0011674 0.04291508 -0.71993332
## 6 -0.671261233 0.5743076 -1.17050260 1.01329858
## logKpre_apap_sul
## 1 -1.50864561
## 2 0.48599584
## 3 -0.05666556
## 4 0.78554977
## 5 -0.54997601
## 6 -1.66492170
## logAPAPCYS_HLM_CL: 8.15992019166072e-05
## logAPAPCYS_Km: 0.000967185735063076
## logAPAPGLU_HLM_CL: 0.00185811908010114
## logAPAPGLU_Km: 0.0181986967501468
## logAPAPSUL_HLM_CL: 0.00194438305778717
## logKa_apap: 0.0052710585929262
## logKpki_apap: 1.14556190846099
## logKpki_apap_cys: 0.627517333300021
## logKpki_apap_glu: 1.72115983885983
## logKpli_apap: 0.686745011371109
## logKpli_apap_cys: 0.741951487277269
## logKpli_apap_glu: 1.09832748267348
## logKpre_apap: 0.954180830717175
## logKpre_apap_cys: 8.93813625877344
## logKpre_apap_glu: 0.277980630583032
## logKpre_apap_sul: 0.22120937811377
In comparison to try 6, I also drop column n before converting the data to a datalist. I use the best fit from try 6 as center for the sampling
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters6 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap", #"F_apap_sul" ,
"Kpre_apap", "Kpki_apap", "Kpli_apap",
"Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
"Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
"Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters6 <- pars[!(names(pars)%in%c(free_parameters6,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters6] <- paste0("exp(log", free_parameters6, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p6_1 <- NULL
for(i in 1:length(p_list)) { p6_1 <<- p6_1 + p_list[[i]]}
pouter6_1 <- log(pars[free_parameters6]) %>% set_names(paste0("log",names(.)))
best_fit <- myfit6 %>% as.parframe() %>% {.[2,]} %>% as.parvec() # der 2. sieht viel besser aus von der dynamik her
pouter6_1[names(best_fit)] <- best_fit
# mypred <- (g*x*p6_1)(seq(0, 48*3600, length.out = 200), pouter6_1, deriv = F)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj6_1 <- normL2(mydatalist, (g*x*p6_1))
# job6_1 <- runbg({myfit <- mstrust(objfun = obj6_1, center = pouter6_1, studyname = "methacetin", cores = 12, fits = 200, iterlim = 200, sd = 2); myfit}, machine = "knecht5", filename = "job6_1", input = global_env_without(c("mypred", "job", "myfit")))
# save(job6_1, file = "job6_1.rda")
# job6_1$check()
# load("job6_1.rda")
# myfit6_1 <- job6_1$get()$knecht5
# save(myfit6_1, file = "myfit6_1.rda")
# job6_1$purge()
load("myfit6_1.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit6_1 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit6_1 %>% as.parframe())+scale_y_log10()
mypred6_1 <- (g*x*p6_1)(mytimes*4, myfit6_1 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot6_1 <- plotCombined(mypred6_1, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot6_1)
# myplot6_1
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 126 2254.796 TRUE 168 -7.160945 -4.394436
## 2 55 2616.813 FALSE 200 -6.451695 -7.616642
## 3 17 2724.360 TRUE 107 5.557814 6.673239
## 4 40 3089.612 TRUE 86 4.825572 9.344937
## 5 49 3140.781 FALSE 200 -7.690501 -8.975971
## 6 152 3143.196 TRUE 57 -7.021284 -1.557577
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -6.718847 -8.311169 -1.954417 -6.524918
## 2 -2.815882 -9.898235 -11.899718 -4.243046
## 3 -4.036551 -9.848632 -6.435657 -7.323644
## 4 -5.141421 -9.827778 -10.177033 -6.097087
## 5 -4.286810 -7.639483 -2.995771 -6.741065
## 6 -7.080733 -10.818246 -6.381171 -5.752037
## logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys logKpki_apap_cys
## 1 0.011772621 1.1690414 0.4379879 0.9132539 -1.9017606
## 2 -0.217629833 1.1177705 -3.9170836 0.4071650 1.7266591
## 3 -0.270321088 -1.2532600 -2.0291797 0.6635614 -1.2993619
## 4 -0.296894373 2.2177932 -0.7600494 -1.0904734 2.5744728
## 5 0.005178141 1.4576398 -2.2549729 -0.7537454 -2.2248627
## 6 0.234766940 0.9271649 1.8868622 -1.9757345 0.2620245
## logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu logKpli_apap_glu
## 1 -3.2571647 -1.58858870 0.4221563 -0.4352780
## 2 0.0398897 -1.46330331 1.4233708 -0.4196139
## 3 -2.7804837 0.05884825 -0.1986642 -1.6136809
## 4 1.5410468 -3.05246767 -2.0834014 -4.2756599
## 5 2.8185400 -1.92875393 0.2224316 5.5323154
## 6 -1.5929422 -3.93035938 -0.3190701 -3.6514846
## logKpre_apap_sul
## 1 -0.6014329
## 2 -0.4671365
## 3 0.2192896
## 4 1.7863395
## 5 -0.4661242
## 6 1.1538600
## logAPAPCYS_HLM_CL: 0.000245756553408302
## logAPAPCYS_Km: 0.141647012033554
## logAPAPGLU_HLM_CL: 0.000776320239293284
## logAPAPGLU_Km: 0.0123458440700416
## logAPAPSUL_HLM_CL: 0.00120793042044707
## logKa_apap: 0.0014664396901935
## logKpki_apap: 3.21890543991051
## logKpki_apap_cys: 0.1493055246381
## logKpki_apap_glu: 1.52524686894739
## logKpli_apap: 1.54958613737682
## logKpli_apap_cys: 0.038497394575441
## logKpli_apap_glu: 0.647084766249551
## logKpre_apap: 1.01184219056339
## logKpre_apap_cys: 2.4924194943794
## logKpre_apap_glu: 0.204213614409576
## logKpre_apap_sul: 0.548025794974965
global_env_without <- function(reg) ls(.GlobalEnv)[!(ls(.GlobalEnv) %>% sapply(. %>% str_detect(reg) %>% any))]
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters7 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters7 <- pars[!(names(pars)%in%c(free_parameters7,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables7 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters7 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters7 <- c(free_parameters7, scale_parameters7)
i <- 2
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters7] <- paste0("exp(log", free_parameters7, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters7, x = scale_parameters7, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters7)
scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any %>% not)] <- "1"
trafo <- c(trafo, scales)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p7 <- NULL
for(i in 1:length(p_list)) { p7 <<- p7 + p_list[[i]]}
g7 <- Y(observables7, x)#, parameters = c(free_parameters7, scale_parameters7))
obj7 <- normL2(mydatalist, (g7*x*p7))
pouter7 <- rep(0, length(getParameters(obj7))) %>% set_names(getParameters(obj7))
pouter7[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# job7 <- runbg({myfit <- mstrust(objfun = obj7, center = pouter7, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht5", filename = "job7", input = global_env_without(c("mypred", "job", "myfit")))
# save(job7, file = "job7.rda")
# job7$check()
# myfit7 <- job7$get()$knecht5
# save(myfit7, file = "myfit7.rda")
# job7$purge()
load("myfit7.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit7 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit7 %>% as.parframe())+scale_y_log10()
mypred7 <- (g7*x*p7)(mytimes*4, myfit7 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot7 <- plotCombined(mypred7, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot7)
# myplot7
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 5 2490.012 TRUE 30 -5.749303 -3.08097475
## 2 64 2626.769 TRUE 60 -4.454476 -2.31516013
## 3 67 2855.601 TRUE 87 -5.293002 -2.39553148
## 4 35 3037.778 FALSE 100 -4.055656 -2.26447288
## 5 98 3347.208 FALSE 100 -5.763790 -4.18377196
## 6 14 3572.777 FALSE 100 -2.353277 -0.00949594
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -5.069446 -10.95827 -6.696778 -7.438818
## 2 -5.574571 -12.31058 -6.987258 -7.603783
## 3 -5.345860 -11.73765 -6.328117 -6.926927
## 4 -5.365076 -11.24999 -7.009670 -7.954771
## 5 -5.881820 -10.51708 -7.158013 -8.062787
## 6 -5.624627 -10.82433 -6.874139 -7.609710
## log_scale_apap_Chan1997_1.4_0 log_scale_apap_glu_Chan1997_1.4_0
## 1 0.7357347 -0.258980454
## 2 0.9075704 0.003338643
## 3 0.4662637 0.365011207
## 4 1.2621657 -0.338130480
## 5 1.0203440 -0.713985018
## 6 0.7625238 0.359423877
## log_scale_apap_sul_Chan1997_1.4_0 log_scale_apap_cys_Chan1997_1.4_0
## 1 -0.03917457 0.6911870
## 2 0.21839269 0.4441861
## 3 -0.26185349 1.4962657
## 4 0.48438243 0.9564866
## 5 0.35501653 0.2993361
## 6 -0.08656187 0.4475230
## log_scale_apap_Chiew2010_5.6_0 log_scale_apap_glu_Chiew2010_5.6_0
## 1 0.28907910 1.4721450
## 2 0.32685873 0.8062451
## 3 0.04167644 0.5089211
## 4 0.28628828 0.1863739
## 5 -0.25747805 -0.5369452
## 6 0.60731079 -1.4552898
## log_scale_apap_sul_Chiew2010_5.6_0 log_scale_apap_Critchley2005_1.4_0
## 1 -0.1620631 -0.3123042
## 2 0.8175662 0.1778291
## 3 -2.1840586 -0.8736844
## 4 0.8380141 0.2651291
## 5 0.2400412 -1.6486526
## 6 0.6459171 0.1341244
## log_scale_apap_glu_Critchley2005_1.4_0
## 1 0.6340540
## 2 -0.8673394
## 3 -0.5829681
## 4 0.1720526
## 5 -0.6356093
## 6 -1.1555520
## log_scale_apap_sul_Critchley2005_1.4_0
## 1 -0.32952950
## 2 -0.58790653
## 3 -0.75363120
## 4 -0.42899409
## 5 0.36572246
## 6 0.06592612
## log_scale_apap_cys_Critchley2005_1.4_0 log_scale_apap_Rawlins1977_0_1
## 1 0.4520727 0.12354755
## 2 1.9508658 0.19603601
## 3 1.1278751 0.01086081
## 4 0.8794284 0.34127665
## 5 -0.4761007 0.23294843
## 6 0.3483191 0.15601077
## log_scale_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 0.27999826 0.23205478
## 2 0.40194965 0.40054980
## 3 0.07005591 0.03746354
## 4 0.69124502 0.72483500
## 5 0.55520067 0.47151149
## 6 0.28847757 0.35586546
## log_scale_apap_Rawlins1977_0.5_0
## 1 -0.32097913
## 2 -0.03674230
## 3 -0.03173343
## 4 -0.78815698
## 5 0.17669530
## 6 0.45105346
## logAPAPCYS_HLM_CL: 1.74133804776165e-05
## logAPAPCYS_Km: 0.00123488444599134
## logAPAPGLU_HLM_CL: 0.0031850000915229
## logAPAPGLU_Km: 0.0459144794614278
## logAPAPSUL_HLM_CL: 0.00628590132812984
## logKa_apap: 0.000587979808605786
## log_scale_apap_Chan1997_1.4_0: 2.08701470180148
## log_scale_apap_Chiew2010_5.6_0: 1.3351973371447
## log_scale_apap_Critchley2005_1.4_0: 0.731758900597009
## log_scale_apap_cys_Chan1997_1.4_0: 1.99608345214162
## log_scale_apap_cys_Critchley2005_1.4_0: 1.57156625159931
## log_scale_apap_glu_Chan1997_1.4_0: 0.771838109420778
## log_scale_apap_glu_Chiew2010_5.6_0: 4.35857412292258
## log_scale_apap_glu_Critchley2005_1.4_0: 1.8852379216948
## log_scale_apap_Rawlins1977_0_1: 1.13150380172748
## log_scale_apap_Rawlins1977_0.5_0: 0.725438390393933
## log_scale_apap_Rawlins1977_1_0: 1.32312750633397
## log_scale_apap_Rawlins1977_2_0: 1.26118881349357
## log_scale_apap_sul_Chan1997_1.4_0: 0.96158283149328
## log_scale_apap_sul_Chiew2010_5.6_0: 0.850387509055339
## log_scale_apap_sul_Critchley2005_1.4_0: 0.719262064522545
Very bad results… Mistake: Don’t remove the rows with sigma=NA! Stupid…
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters8 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters8 <- pars[!(names(pars)%in%c(free_parameters8,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables8 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters8 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters8 <- c(free_parameters8, scale_parameters8)
error_model8 <- c(apap = "srel_apap*apap^2 +s0_apap",
apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
error_parameters8 <- setdiff(getSymbols(error_model8), names(error_model8)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters8] <- paste0("exp(log", free_parameters8, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters8, x = scale_parameters8, y = .)}
scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters8, x = error_parameters8, y = .)}
errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p8 <- NULL
for(i in 1:length(p_list)) { p8 <<- p8 + p_list[[i]]}
g8 <- Y(observables8, x)#, parameters = c(free_parameters8, scale_parameters8))
err8 <- Y(error_model8, g8)
obj8 <- normL2(mydatalist, (g8*x*p8), errmodel = err8)
pouter8 <- rep(0, length(getParameters(obj8))) %>% set_names(getParameters(obj8))
pouter8[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# obj8(pouter8)
# job8 <- runbg({myfit <- mstrust(objfun = obj8, center = pouter8, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht5", filename = "job8", input = global_env_without(c("mypred", "job", "myfit")))
# save(job8, file = "job8.rda")
# job8$purge()
# myfit8 <- job8$get()$knecht5
# save(myfit8, file = "myfit8.rda")
# job8$purge()
load("myfit8.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit8 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit8 %>% as.parframe())#+scale_y_log10()
mypred8 <- (g8*x*p8)(mytimes*4, myfit8 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot8 <- plotCombined(mypred8, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot8)
# myplot8
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 46 -486.5434 TRUE 68 0.8448040 3.0690138
## 2 95 -477.9974 TRUE 86 2.3221006 4.5447823
## 3 3 -458.5112 TRUE 85 -2.8788000 -0.5052364
## 4 91 -451.1060 TRUE 84 2.0359641 4.2579075
## 5 50 -437.6416 FALSE 100 2.1876617 4.4051112
## 6 44 -426.7063 TRUE 82 -0.9089178 1.3108703
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -10.788256 -13.56917 -3.126300 -6.295179
## 2 -12.281700 -11.21214 -1.128100 -6.280291
## 3 -6.651499 -14.82540 -3.552027 -6.325155
## 4 -12.473814 -13.54999 -2.420064 -6.289021
## 5 -11.711911 -13.73429 -2.595443 -2.278541
## 6 -10.442107 -14.43848 -3.406845 -2.307167
## log_scale_apap_Chan1997_1.4_0 log_scale_apap_glu_Chan1997_1.4_0
## 1 -0.34147420 0.3524540
## 2 0.87395974 -0.2202193
## 3 0.28419894 -1.2117128
## 4 -0.23704451 -0.5141464
## 5 -0.00256986 0.7640506
## 6 -1.39510256 -0.8312082
## log_scale_apap_sul_Chan1997_1.4_0 log_scale_apap_cys_Chan1997_1.4_0
## 1 -0.09054225 0.4909674
## 2 -2.05176756 0.1713548
## 3 1.59152206 0.9567482
## 4 -0.64296093 0.6822799
## 5 -0.28381384 1.2033634
## 6 -0.42447087 -0.9383802
## log_srel_apap_Chan1997_1.4_0 log_s0_apap_Chan1997_1.4_0
## 1 0.2569589 -1.9014794
## 2 -0.6194964 1.4224062
## 3 0.4884257 0.4620627
## 4 0.4431686 0.3437440
## 5 -1.0774869 -0.7111688
## 6 -0.1059288 -0.4287323
## log_srel_apap_glu_Chan1997_1.4_0 log_s0_apap_glu_Chan1997_1.4_0
## 1 -0.9569274 -0.1547545
## 2 0.6830946 -0.9358321
## 3 0.1252113 1.4079458
## 4 -1.0505252 -1.1969937
## 5 1.3007331 -1.5656455
## 6 2.1081213 -1.7762023
## log_srel_apap_sul_Chan1997_1.4_0 log_s0_apap_sul_Chan1997_1.4_0
## 1 -2.5583855 -1.1857697
## 2 0.2056214 -1.4048882
## 3 -0.5206815 -3.0952762
## 4 1.7018681 -2.4248608
## 5 -1.9290153 -0.8441127
## 6 0.6754411 0.6658790
## log_srel_apap_cys_Chan1997_1.4_0 log_s0_apap_cys_Chan1997_1.4_0
## 1 -0.80413329 -1.76384101
## 2 -1.03467092 -2.63912882
## 3 0.12948514 -0.06831229
## 4 -1.23494485 -0.35465510
## 5 -1.45640160 -2.21684171
## 6 -0.08965077 -2.18967728
## log_scale_apap_Chiew2010_5.6_0 log_scale_apap_glu_Chiew2010_5.6_0
## 1 -0.73225919 0.29875836
## 2 0.07746059 -0.75638815
## 3 -0.41557105 0.66860572
## 4 -0.34446349 -0.09935669
## 5 -0.06790696 0.17262465
## 6 -0.36695303 0.98236576
## log_scale_apap_sul_Chiew2010_5.6_0 log_srel_apap_Chiew2010_5.6_0
## 1 -0.2751353 0.2001658
## 2 -0.8853731 -1.1763805
## 3 -1.4579817 1.4459862
## 4 1.9013839 -1.1825723
## 5 0.8957957 0.4329237
## 6 -0.5285034 -1.6273736
## log_s0_apap_Chiew2010_5.6_0 log_srel_apap_glu_Chiew2010_5.6_0
## 1 -1.3083801 0.2537679
## 2 -0.9910735 0.4781154
## 3 -1.3495770 -0.2456441
## 4 -0.4922068 -0.4925589
## 5 1.2760837 1.4551869
## 6 -1.4419611 -2.2427514
## log_s0_apap_glu_Chiew2010_5.6_0 log_srel_apap_sul_Chiew2010_5.6_0
## 1 -0.9917074 0.30899346
## 2 -0.1322426 0.20086298
## 3 -2.4500674 1.13602514
## 4 -2.1513778 0.97349583
## 5 0.3285409 0.09686059
## 6 -1.1718316 -0.32388955
## log_s0_apap_sul_Chiew2010_5.6_0 log_scale_apap_Critchley2005_1.4_0
## 1 -0.2317920 0.75661973
## 2 -0.1080672 0.40330247
## 3 -1.2538891 -1.84876130
## 4 0.5970107 0.09984706
## 5 0.3377952 -1.08319140
## 6 -0.9579787 0.15268520
## log_scale_apap_glu_Critchley2005_1.4_0
## 1 -1.6075093
## 2 -0.4203937
## 3 -0.4089714
## 4 0.6322839
## 5 0.5396460
## 6 0.6909803
## log_scale_apap_sul_Critchley2005_1.4_0
## 1 1.1200586
## 2 1.3772116
## 3 -0.8219418
## 4 -0.1005233
## 5 -1.1104974
## 6 -0.4086092
## log_scale_apap_cys_Critchley2005_1.4_0 log_srel_apap_Critchley2005_1.4_0
## 1 0.7914627 0.5229484
## 2 -0.2206032 1.7093085
## 3 -0.7028494 -0.4162756
## 4 -0.3442296 -1.5224635
## 5 -0.7764242 -1.5065604
## 6 -1.2219040 -0.4096500
## log_s0_apap_Critchley2005_1.4_0 log_srel_apap_glu_Critchley2005_1.4_0
## 1 2.8885995 -0.6958732
## 2 0.2769623 1.6113013
## 3 1.9668600 0.8816571
## 4 0.4696681 0.6190646
## 5 -0.1120432 0.9772918
## 6 2.0007139 -1.2027400
## log_s0_apap_glu_Critchley2005_1.4_0
## 1 -1.462525446
## 2 -0.876828265
## 3 -2.072184801
## 4 -0.374247260
## 5 -0.093214581
## 6 0.005710159
## log_srel_apap_sul_Critchley2005_1.4_0
## 1 0.09474047
## 2 1.05831194
## 3 1.37006018
## 4 0.49068959
## 5 -0.15389290
## 6 -0.21423115
## log_s0_apap_sul_Critchley2005_1.4_0
## 1 -1.5436022
## 2 -1.9787239
## 3 1.1611746
## 4 -0.5722336
## 5 -0.7604901
## 6 -0.5263634
## log_srel_apap_cys_Critchley2005_1.4_0
## 1 0.6346692
## 2 -0.9788842
## 3 -1.7831523
## 4 0.6156580
## 5 -1.3146436
## 6 0.4575737
## log_s0_apap_cys_Critchley2005_1.4_0 log_scale_apap_Rawlins1977_0_1
## 1 0.4977260 0.021621889
## 2 -0.6198171 0.021845011
## 3 -1.8375393 0.008699087
## 4 -1.2241609 0.022133087
## 5 -1.1091981 0.026012429
## 6 -0.1110733 0.024257471
## log_srel_apap_Rawlins1977_0_1 log_s0_apap_Rawlins1977_0_1
## 1 2.454559 -8.977730
## 2 2.454535 -8.977814
## 3 2.483257 -8.996088
## 4 2.453265 -8.974551
## 5 2.437694 -8.941497
## 6 2.442176 -8.947284
## log_scale_apap_Rawlins1977_1_0 log_srel_apap_Rawlins1977_1_0
## 1 0.1678273 2.582696
## 2 0.1692129 2.581302
## 3 0.1388991 2.595842
## 4 0.1692214 2.580949
## 5 0.2307545 2.589881
## 6 0.2203842 2.591367
## log_s0_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 -14.66563 0.4137086
## 2 -16.05141 0.4157914
## 3 -15.53460 0.3500338
## 4 -15.71205 0.4157741
## 5 -16.02509 0.4775457
## 6 -15.33039 0.4570333
## log_srel_apap_Rawlins1977_2_0 log_s0_apap_Rawlins1977_2_0
## 1 2.341414 -15.33500
## 2 2.341451 -16.63546
## 3 2.347062 -16.04883
## 4 2.341416 -16.30266
## 5 2.346485 -16.51624
## 6 2.347016 -15.82922
## log_scale_apap_Rawlins1977_0.5_0 log_srel_apap_Rawlins1977_0.5_0
## 1 -1.5900253 1.9345084
## 2 -0.2137423 0.2236760
## 3 0.8627395 0.6101951
## 4 1.1149818 -0.5099961
## 5 -0.4089437 1.0647446
## 6 -0.1155511 0.7124745
## log_s0_apap_Rawlins1977_0.5_0
## 1 -2.5199219
## 2 -0.2345078
## 3 -1.2819850
## 4 1.5677611
## 5 -0.7884726
## 6 -0.3613277
## logAPAPCYS_HLM_CL: 1.2793391679861e-06
## logAPAPCYS_Km: 0.0438798488903149
## logAPAPGLU_HLM_CL: 2.32752164063581
## logAPAPGLU_Km: 21.520668349057
## logAPAPSUL_HLM_CL: 2.06404786189464e-05
## logKa_apap: 0.00184517865281865
## log_s0_apap_Chan1997_1.4_0: 0.149347509305303
## log_s0_apap_Chiew2010_5.6_0: 0.270257502535883
## log_s0_apap_Critchley2005_1.4_0: 17.9681281235639
## log_s0_apap_cys_Chan1997_1.4_0: 0.171385306079452
## log_s0_apap_cys_Critchley2005_1.4_0: 1.64497626822607
## log_s0_apap_glu_Chan1997_1.4_0: 0.856625476079955
## log_s0_apap_glu_Chiew2010_5.6_0: 0.370942798641939
## log_s0_apap_glu_Critchley2005_1.4_0: 0.231650514577204
## log_s0_apap_Rawlins1977_0_1: 0.000126188936947135
## log_s0_apap_Rawlins1977_0.5_0: 0.0804658932475238
## log_s0_apap_Rawlins1977_1_0: 4.27362459933834e-07
## log_s0_apap_Rawlins1977_2_0: 2.18822833245219e-07
## log_s0_apap_sul_Chan1997_1.4_0: 0.30551094502221
## log_s0_apap_sul_Chiew2010_5.6_0: 0.793111039544634
## log_s0_apap_sul_Critchley2005_1.4_0: 0.21361025170362
## log_scale_apap_Chan1997_1.4_0: 0.710721805251837
## log_scale_apap_Chiew2010_5.6_0: 0.480821497303708
## log_scale_apap_Critchley2005_1.4_0: 2.13106047495752
## log_scale_apap_cys_Chan1997_1.4_0: 1.63389603961365
## log_scale_apap_cys_Critchley2005_1.4_0: 2.20662164565056
## log_scale_apap_glu_Chan1997_1.4_0: 1.42255421536611
## log_scale_apap_glu_Chiew2010_5.6_0: 1.34818380327456
## log_scale_apap_glu_Critchley2005_1.4_0: 0.200386104365974
## log_scale_apap_Rawlins1977_0_1: 1.02185733567333
## log_scale_apap_Rawlins1977_0.5_0: 0.203920455093923
## log_scale_apap_Rawlins1977_1_0: 1.18273235599321
## log_scale_apap_Rawlins1977_2_0: 1.51241638200544
## log_scale_apap_sul_Chan1997_1.4_0: 0.913435744877439
## log_scale_apap_sul_Chiew2010_5.6_0: 0.759469382839319
## log_scale_apap_sul_Critchley2005_1.4_0: 3.06503365611115
## log_srel_apap_Chan1997_1.4_0: 1.29299192947592
## log_srel_apap_Chiew2010_5.6_0: 1.22160522476788
## log_srel_apap_Critchley2005_1.4_0: 1.68699428030717
## log_srel_apap_cys_Chan1997_1.4_0: 0.447475592070443
## log_srel_apap_cys_Critchley2005_1.4_0: 1.88639794692755
## log_srel_apap_glu_Chan1997_1.4_0: 0.384071160069748
## log_srel_apap_glu_Chiew2010_5.6_0: 1.28887257124136
## log_srel_apap_glu_Critchley2005_1.4_0: 0.498638857365226
## log_srel_apap_Rawlins1977_0_1: 11.6412981287543
## log_srel_apap_Rawlins1977_0.5_0: 6.92064127305638
## log_srel_apap_Rawlins1977_1_0: 13.2327664425859
## log_srel_apap_Rawlins1977_2_0: 10.3959247377831
## log_srel_apap_sul_Chan1997_1.4_0: 0.0774296530505137
## log_srel_apap_sul_Chiew2010_5.6_0: 1.36205346050779
## log_srel_apap_sul_Critchley2005_1.4_0: 1.09937349942295
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters9 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters9 <- pars[!(names(pars)%in%c(free_parameters9,names(f)[1]))] %>% names
mydatalist <- data %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables9 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters9 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters9 <- c(free_parameters9, scale_parameters9)
error_model9 <- c(apap = "srel_apap*apap^2 +s0_apap",
apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
error_parameters9 <- setdiff(getSymbols(error_model9), names(error_model9)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters9] <- paste0("exp(log", free_parameters9, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters9, x = scale_parameters9, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters9)
scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters9, x = error_parameters9, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters9)
errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p9 <- NULL
for(i in 1:length(p_list)) { p9 <<- p9 + p_list[[i]]}
g9 <- Y(observables9, x)#, parameters = c(free_parameters9, scale_parameters9))
err9 <- Y(error_model9, g9)
obj9 <- normL2(mydatalist, (g9*x*p9), errmodel = err9)
pouter9 <- rep(0, length(getParameters(obj9))) %>% set_names(getParameters(obj9))
pouter9[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# obj9(pouter9)
# job9 <- runbg({myfit <- mstrust(objfun = obj9, center = pouter9, studyname = "methacetin", cores = 16, fits = 200); myfit}, machine = "ruprecht2", filename = "job9", input = global_env_without(c("mypred", "job", "myfit")))
# save(job9, file = "job9.rda")
# job9$check()
#
# load("job9.rda")
# myfit9 <- job9$get()$ruprecht2
# save(myfit9, file = "myfit9.rda")
# job9$purge()
load("myfit9.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit9 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit9 %>% as.parframe())#+scale_y_log10()
mypred9 <- (g9*x*p9)(mytimes*4, myfit9 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot9 <- plotCombined(mypred9, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot9)
# myplot9
# obj9(myfit9 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.})
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 1 -2477.716 FALSE 100 -5.099331 7.999288
## 2 35 -2477.716 FALSE 100 -5.032580 7.301881
## 3 49 -2477.715 FALSE 100 -5.087414 7.628175
## 4 179 -2477.715 FALSE 100 -5.087632 7.116519
## 5 64 -2477.715 TRUE 99 -5.073408 7.132834
## 6 127 -2477.715 FALSE 100 -5.094321 6.967860
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -4.797299 -17.83559 -6.250228 -7.073747
## 2 -4.797371 -17.08140 -6.250285 -7.073740
## 3 -4.797351 -17.45008 -6.250078 -7.073678
## 4 -4.797346 -16.93769 -6.250345 -7.073796
## 5 -4.797347 -16.94634 -6.250387 -7.073790
## 6 -4.797427 -16.79246 -6.250201 -7.073706
## log_scale_apap_Chan1997_14_0 log_scale_apap_glu_Chan1997_14_0
## 1 0.3201883 10.240255
## 2 0.3201823 9.476110
## 3 0.3201498 9.857194
## 4 0.3202139 9.345829
## 5 0.3202098 9.347900
## 6 0.3201586 9.203812
## log_scale_apap_sul_Chan1997_14_0 log_scale_apap_cys_Chan1997_14_0
## 1 -0.5703412 7.541992
## 2 -0.5702552 6.787785
## 3 -0.5703223 7.156501
## 4 -0.5702481 6.644072
## 5 -0.5702448 6.652709
## 6 -0.5702149 6.498853
## log_srel_apap_Chan1997_14_0 log_s0_apap_Chan1997_14_0
## 1 -1.066589 -8.553130
## 2 -1.066910 -8.552872
## 3 -1.064453 -8.553918
## 4 -1.069107 -8.552475
## 5 -1.070049 -8.552360
## 6 -1.066075 -8.552996
## log_srel_apap_glu_Chan1997_14_0 log_s0_apap_glu_Chan1997_14_0
## 1 3.937521 -7.368639
## 2 3.937853 -7.368713
## 3 3.937651 -7.368663
## 4 3.937780 -7.368695
## 5 3.937766 -7.368738
## 6 3.937715 -7.368682
## log_srel_apap_sul_Chan1997_14_0 log_s0_apap_sul_Chan1997_14_0
## 1 4.178282 -9.337702
## 2 4.178646 -9.337195
## 3 4.178415 -9.337887
## 4 4.178152 -9.337426
## 5 4.178408 -9.337071
## 6 4.178751 -9.337377
## log_srel_apap_cys_Chan1997_14_0 log_s0_apap_cys_Chan1997_14_0
## 1 -4.014309 -9.955108
## 2 -3.320616 -9.955241
## 3 -3.635824 -9.955119
## 4 -3.161390 -9.955246
## 5 -3.162397 -9.955206
## 6 -2.983036 -9.955250
## log_scale_apap_Chiew2010_56_0 log_scale_apap_glu_Chiew2010_56_0
## 1 0.1542391 10.561511
## 2 0.1542455 9.797452
## 3 0.1541823 10.178496
## 4 0.1542872 9.667135
## 5 0.1542769 9.669024
## 6 0.1542120 9.525162
## log_scale_apap_sul_Chiew2010_56_0 log_srel_apap_Chiew2010_56_0
## 1 -0.7426937 0.3778347
## 2 -0.7426260 0.3779957
## 3 -0.7426694 0.3771786
## 4 -0.7426272 0.3780723
## 5 -0.7426270 0.3777571
## 6 -0.7425874 0.3779534
## log_s0_apap_Chiew2010_56_0 log_srel_apap_glu_Chiew2010_56_0
## 1 -7.069136 2.369965
## 2 -7.067144 2.370166
## 3 -7.067803 2.370305
## 4 -7.066808 2.371076
## 5 -7.066785 2.371267
## 6 -7.067703 2.370746
## log_s0_apap_glu_Chiew2010_56_0 log_srel_apap_sul_Chiew2010_56_0
## 1 -5.863964 1.413941
## 2 -5.863894 1.413643
## 3 -5.864205 1.414305
## 4 -5.864437 1.413383
## 5 -5.865029 1.413175
## 6 -5.864204 1.413946
## log_s0_apap_sul_Chiew2010_56_0 log_scale_apap_Critchley2005_14_0
## 1 -6.186949 0.2345061
## 2 -6.187438 0.2344554
## 3 -6.187013 0.2343445
## 4 -6.187419 0.2344923
## 5 -6.187215 0.2345560
## 6 -6.187368 0.2344314
## log_scale_apap_glu_Critchley2005_14_0
## 1 10.314313
## 2 9.550148
## 3 9.931238
## 4 9.419870
## 5 9.421957
## 6 9.277840
## log_scale_apap_sul_Critchley2005_14_0
## 1 -1.148680
## 2 -1.148614
## 3 -1.148658
## 4 -1.148615
## 5 -1.148612
## 6 -1.148576
## log_scale_apap_cys_Critchley2005_14_0 log_srel_apap_Critchley2005_14_0
## 1 7.232789 2.288037
## 2 6.478583 2.284297
## 3 6.847288 2.279003
## 4 6.334869 2.284317
## 5 6.343503 2.285165
## 6 6.189646 2.284663
## log_s0_apap_Critchley2005_14_0 log_srel_apap_glu_Critchley2005_14_0
## 1 -8.080480 -5.922499
## 2 -8.074298 -5.172423
## 3 -8.068500 -5.541084
## 4 -8.073987 -5.031696
## 5 -8.075026 -5.034967
## 6 -8.075041 -4.887071
## log_s0_apap_glu_Critchley2005_14_0 log_srel_apap_sul_Critchley2005_14_0
## 1 -6.304482 4.463234
## 2 -6.304953 4.464450
## 3 -6.304470 4.464503
## 4 -6.304871 4.464421
## 5 -6.304808 4.464391
## 6 -6.304898 4.464320
## log_s0_apap_sul_Critchley2005_14_0 log_srel_apap_cys_Critchley2005_14_0
## 1 -7.946701 -2.811800
## 2 -7.945526 -2.085599
## 3 -7.945487 -2.428260
## 4 -7.945550 -1.991045
## 5 -7.945590 -1.941678
## 6 -7.945679 -1.917429
## log_s0_apap_cys_Critchley2005_14_0 log_scale_apap_Rawlins1977_0_1
## 1 -10.12659 -0.1015749
## 2 -10.12657 -0.1015765
## 3 -10.12661 -0.1016133
## 4 -10.12658 -0.1015472
## 5 -10.12661 -0.1015634
## 6 -10.12657 -0.1015950
## log_srel_apap_Rawlins1977_0_1 log_s0_apap_Rawlins1977_0_1
## 1 2.881313 -10.45464
## 2 2.880629 -10.44504
## 3 2.880721 -10.44162
## 4 2.880528 -10.44397
## 5 2.879162 -10.41924
## 6 2.880092 -10.43567
## log_scale_apap_Rawlins1977_1_0 log_srel_apap_Rawlins1977_1_0
## 1 -0.03506085 3.061212
## 2 -0.03507190 3.062111
## 3 -0.03511765 3.062435
## 4 -0.03503499 3.062118
## 5 -0.03503313 3.062624
## 6 -0.03510720 3.062807
## log_s0_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 -8.476342 0.03209758
## 2 -8.477932 0.03208771
## 3 -8.478216 0.03205411
## 4 -8.478040 0.03212816
## 5 -8.478879 0.03214222
## 6 -8.479337 0.03205324
## log_srel_apap_Rawlins1977_2_0 log_s0_apap_Rawlins1977_2_0
## 1 2.456148 -8.369891
## 2 2.457079 -8.371340
## 3 2.457015 -8.370964
## 4 2.457062 -8.371123
## 5 2.457133 -8.372795
## 6 2.456935 -8.370860
## log_scale_apap_Rawlins1977_05_0 log_srel_apap_Rawlins1977_05_0
## 1 -0.2098564 -7.088126
## 2 -0.2098571 -6.338162
## 3 -0.2098887 -6.706408
## 4 -0.2098285 -6.195749
## 5 -0.2098314 -6.201175
## 6 -0.2098756 -6.053133
## log_s0_apap_Rawlins1977_05_0
## 1 -8.283322
## 2 -8.284305
## 3 -8.283902
## 4 -8.284432
## 5 -8.284290
## 6 -8.284218
## logAPAPCYS_HLM_CL: 1.795146994404e-08
## logAPAPCYS_Km: 0.00193001355192211
## logAPAPGLU_HLM_CL: 0.00610082435223236
## logAPAPGLU_Km: 2978.83487681056
## logAPAPSUL_HLM_CL: 0.00825200823079057
## logKa_apap: 0.000847053225794413
## log_s0_apap_Chan1997_14_0: 0.00019294018495771
## log_s0_apap_Chiew2010_56_0: 0.00085096838208337
## log_s0_apap_Critchley2005_14_0: 0.000309522401143699
## log_s0_apap_cys_Chan1997_14_0: 4.74844786235421e-05
## log_s0_apap_cys_Critchley2005_14_0: 4.00018130602822e-05
## log_s0_apap_glu_Chan1997_14_0: 0.000630726272006818
## log_s0_apap_glu_Chiew2010_56_0: 0.00283996453762311
## log_s0_apap_glu_Critchley2005_14_0: 0.00182809286937598
## log_s0_apap_Rawlins1977_0_1: 2.8814215684246e-05
## log_s0_apap_Rawlins1977_05_0: 0.000252696390036355
## log_s0_apap_Rawlins1977_1_0: 0.000208339338154456
## log_s0_apap_Rawlins1977_2_0: 0.000231740837081231
## log_s0_apap_sul_Chan1997_14_0: 8.80414871966006e-05
## log_s0_apap_sul_Chiew2010_56_0: 0.00205609044641647
## log_s0_apap_sul_Critchley2005_14_0: 0.000353827456947821
## log_scale_apap_Chan1997_14_0: 1.37738706359055
## log_scale_apap_Chiew2010_56_0: 1.16676983632714
## log_scale_apap_Critchley2005_14_0: 1.26428417568878
## log_scale_apap_cys_Chan1997_14_0: 1885.5826917627
## log_scale_apap_cys_Critchley2005_14_0: 1384.0773914102
## log_scale_apap_glu_Chan1997_14_0: 28008.2659198081
## log_scale_apap_glu_Chiew2010_56_0: 38619.4394146001
## log_scale_apap_glu_Critchley2005_14_0: 30161.2433788063
## log_scale_apap_Rawlins1977_0_1: 0.903413505390761
## log_scale_apap_Rawlins1977_05_0: 0.81070064561698
## log_scale_apap_Rawlins1977_1_0: 0.965546657143368
## log_scale_apap_Rawlins1977_2_0: 1.03261825838704
## log_scale_apap_sul_Chan1997_14_0: 0.565332533789827
## log_scale_apap_sul_Chiew2010_56_0: 0.475830435424811
## log_scale_apap_sul_Critchley2005_14_0: 0.317054858672663
## log_srel_apap_Chan1997_14_0: 0.344180473770218
## log_srel_apap_Chiew2010_56_0: 1.45912166820889
## log_srel_apap_Critchley2005_14_0: 9.85557305660906
## log_srel_apap_cys_Chan1997_14_0: 0.0180554247494083
## log_srel_apap_cys_Critchley2005_14_0: 0.0600966934703123
## log_srel_apap_glu_Chan1997_14_0: 51.2912931730463
## log_srel_apap_glu_Chiew2010_56_0: 10.6970164613051
## log_srel_apap_glu_Critchley2005_14_0: 0.00267849706185607
## log_srel_apap_Rawlins1977_0_1: 17.8376851856982
## log_srel_apap_Rawlins1977_05_0: 0.000834960350013315
## log_srel_apap_Rawlins1977_1_0: 21.3534269462142
## log_srel_apap_Rawlins1977_2_0: 11.6598121177031
## log_srel_apap_sul_Chan1997_14_0: 65.2536831006089
## log_srel_apap_sul_Chiew2010_56_0: 4.11213078704973
## log_srel_apap_sul_Critchley2005_14_0: 86.7676957811533
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters10 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters10 <- pars[!(names(pars)%in%c(free_parameters10,names(f)[1]))] %>% names
mydatalist <- data %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables10 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters10 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters10 <- c(free_parameters10, scale_parameters10)
error_model10 <- c(apap = "srel_apap*apap^2 +s0_apap",
apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
error_parameters10 <- setdiff(getSymbols(error_model10), names(error_model10)) %>% set_names(.,.)
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters10] <- paste0("exp(log", free_parameters10, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters10, x = scale_parameters10, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters10)
scales[1:length(scales)] <- "1"
errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters10, x = error_parameters10, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters10)
errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p10 <- NULL
for(i in 1:length(p_list)) { p10 <<- p10 + p_list[[i]]}
g10 <- Y(observables10, x)#, parameters = c(free_parameters10, scale_parameters10))
err10 <- Y(error_model10, g10)
obj10 <- normL2(mydatalist, (g10*x*p10), errmodel = err10)
pouter10 <- rep(0, length(getParameters(obj10))) %>% set_names(getParameters(obj10))
pouter10[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# job10 <- runbg({myfit <- mstrust(objfun = obj10, center = pouter10, studyname = "methacetin", cores = 16, fits = 200); myfit}, machine = "ruprecht1", filename = "job10", input = global_env_without(c("mypred", "job", "myfit")))
# save(job10, file = "job10.rda")
# job10$check()
# load("job10.rda")
# myfit10 <- job10$get()$ruprecht1
# save(myfit10, file = "myfit10.rda")
# job10$purge()
load("myfit10.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit10 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit10 %>% as.parframe())#+scale_y_log10()
mypred10 <- (g10*x*p10)(mytimes*4, myfit10 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot10 <- plotCombined(mypred10, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot10)
# myplot10
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 72 -2356.519 TRUE 99 -5.984038 -3.407580
## 2 23 -2356.519 TRUE 100 -5.984055 -3.407593
## 3 90 -2356.519 TRUE 98 -5.984039 -3.407588
## 4 144 -2356.519 FALSE 100 -5.984040 -3.407582
## 5 63 -2356.519 TRUE 78 -5.983866 -3.407695
## 6 135 -2356.519 FALSE 100 -5.984081 -3.407616
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -5.782424 -10.57408 -6.696947 -7.013798
## 2 -5.782408 -10.57408 -6.696946 -7.013783
## 3 -5.782429 -10.57408 -6.696951 -7.013733
## 4 -5.782424 -10.57408 -6.696946 -7.013786
## 5 -5.782736 -10.57412 -6.697111 -7.013753
## 6 -5.782389 -10.57409 -6.696947 -7.013793
## log_srel_apap_Chan1997_14_0 log_s0_apap_Chan1997_14_0
## 1 3.289950 -7.052483
## 2 3.289934 -7.052496
## 3 3.289921 -7.052506
## 4 3.289937 -7.052506
## 5 3.289810 -7.052222
## 6 3.289991 -7.052427
## log_srel_apap_glu_Chan1997_14_0 log_s0_apap_glu_Chan1997_14_0
## 1 4.404972 -7.517436
## 2 4.405154 -7.517363
## 3 4.405028 -7.517374
## 4 4.405060 -7.517341
## 5 4.403288 -7.516821
## 6 4.405348 -7.517387
## log_srel_apap_sul_Chan1997_14_0 log_s0_apap_sul_Chan1997_14_0
## 1 4.959255 -8.122914
## 2 4.958957 -8.123157
## 3 4.959125 -8.123015
## 4 4.959096 -8.123068
## 5 4.960016 -8.122510
## 6 4.958673 -8.123290
## log_srel_apap_cys_Chan1997_14_0 log_s0_apap_cys_Chan1997_14_0
## 1 -6.968630 -9.728692
## 2 -3.817880 -9.728688
## 3 -5.081099 -9.728682
## 4 -3.199117 -9.728689
## 5 -2.727214 -9.728656
## 6 -2.482580 -9.728732
## log_srel_apap_Chiew2010_56_0 log_s0_apap_Chiew2010_56_0
## 1 0.1017979 -7.141670
## 2 0.1020375 -7.141873
## 3 0.1016049 -7.141541
## 4 0.1018871 -7.141784
## 5 0.1016913 -7.141588
## 6 0.1019149 -7.141807
## log_srel_apap_glu_Chiew2010_56_0 log_s0_apap_glu_Chiew2010_56_0
## 1 5.497352 -6.479678
## 2 5.497414 -6.479689
## 3 5.497194 -6.479812
## 4 5.497293 -6.479720
## 5 5.496531 -6.479613
## 6 5.498219 -6.479841
## log_srel_apap_sul_Chiew2010_56_0 log_s0_apap_sul_Chiew2010_56_0
## 1 2.464006 -6.105482
## 2 2.463693 -6.105647
## 3 2.463902 -6.105591
## 4 2.463910 -6.105581
## 5 2.468499 -6.106767
## 6 2.463066 -6.105751
## log_srel_apap_Critchley2005_14_0 log_s0_apap_Critchley2005_14_0
## 1 2.350997 -6.825250
## 2 2.350973 -6.825248
## 3 2.350992 -6.825255
## 4 2.350957 -6.825260
## 5 2.350889 -6.824922
## 6 2.350855 -6.825138
## log_srel_apap_glu_Critchley2005_14_0 log_s0_apap_glu_Critchley2005_14_0
## 1 3.602125 -6.967212
## 2 3.602291 -6.967185
## 3 3.602170 -6.967135
## 4 3.602175 -6.967151
## 5 3.601192 -6.966273
## 6 3.602554 -6.967230
## log_srel_apap_sul_Critchley2005_14_0 log_s0_apap_sul_Critchley2005_14_0
## 1 4.017564 -7.716869
## 2 4.017535 -7.716892
## 3 4.017606 -7.716842
## 4 4.017561 -7.716888
## 5 4.017116 -7.717051
## 6 4.017688 -7.716863
## log_srel_apap_cys_Critchley2005_14_0 log_s0_apap_cys_Critchley2005_14_0
## 1 -8.051930 -10.10211
## 2 -4.902004 -10.10212
## 3 -6.164908 -10.10211
## 4 -4.283094 -10.10212
## 5 -3.755340 -10.10221
## 6 -3.566113 -10.10215
## log_srel_apap_Rawlins1977_0_1 log_s0_apap_Rawlins1977_0_1
## 1 2.541371 -8.725966
## 2 2.541383 -8.725981
## 3 2.541391 -8.725948
## 4 2.541383 -8.725955
## 5 2.541252 -8.725793
## 6 2.541297 -8.725894
## log_srel_apap_Rawlins1977_1_0 log_s0_apap_Rawlins1977_1_0
## 1 2.772759 -8.547667
## 2 2.772043 -8.546641
## 3 2.770963 -8.545399
## 4 2.772593 -8.547197
## 5 2.772166 -8.547219
## 6 2.772761 -8.547107
## log_srel_apap_Rawlins1977_2_0 log_s0_apap_Rawlins1977_2_0
## 1 2.485538 -8.345503
## 2 2.485019 -8.343667
## 3 2.485324 -8.344722
## 4 2.485179 -8.344683
## 5 2.485390 -8.344672
## 6 2.485404 -8.344505
## log_srel_apap_Rawlins1977_05_0 log_s0_apap_Rawlins1977_05_0
## 1 3.075085 -8.175010
## 2 3.075391 -8.175087
## 3 3.075176 -8.175017
## 4 3.074992 -8.174969
## 5 3.074871 -8.175119
## 6 3.074939 -8.175055
## logAPAPCYS_HLM_CL: 2.5570243395892e-05
## logAPAPCYS_Km: 0.00123467580981531
## logAPAPGLU_HLM_CL: 0.00251863572465911
## logAPAPGLU_Km: 0.033121273391679
## logAPAPSUL_HLM_CL: 0.00308123765707667
## logKa_apap: 0.000899386426280096
## log_s0_apap_Chan1997_14_0: 0.000865258189660525
## log_s0_apap_Chiew2010_56_0: 0.000791429248666158
## log_s0_apap_Critchley2005_14_0: 0.00108600472415769
## log_s0_apap_cys_Chan1997_14_0: 5.9550106934986e-05
## log_s0_apap_cys_Critchley2005_14_0: 4.09928708157083e-05
## log_s0_apap_glu_Chan1997_14_0: 0.000543524495738245
## log_s0_apap_glu_Chiew2010_56_0: 0.00153430475174298
## log_s0_apap_glu_Critchley2005_14_0: 0.000942275918358174
## log_s0_apap_Rawlins1977_0_1: 0.000162315884677548
## log_s0_apap_Rawlins1977_05_0: 0.000281603560301461
## log_s0_apap_Rawlins1977_1_0: 0.000193997116676406
## log_s0_apap_Rawlins1977_2_0: 0.000237462036719598
## log_s0_apap_sul_Chan1997_14_0: 0.000296662808440626
## log_s0_apap_sul_Chiew2010_56_0: 0.00223060565323054
## log_s0_apap_sul_Critchley2005_14_0: 0.000445252327814248
## log_srel_apap_Chan1997_14_0: 26.8415245012767
## log_srel_apap_Chiew2010_56_0: 1.10715971018415
## log_srel_apap_Critchley2005_14_0: 10.4960243104907
## log_srel_apap_cys_Chan1997_14_0: 0.00094094138090811
## log_srel_apap_cys_Critchley2005_14_0: 0.000318486578636513
## log_srel_apap_glu_Chan1997_14_0: 81.8568173051488
## log_srel_apap_glu_Chiew2010_56_0: 244.044756600744
## log_srel_apap_glu_Critchley2005_14_0: 36.6760946617877
## log_srel_apap_Rawlins1977_0_1: 12.6970698875144
## log_srel_apap_Rawlins1977_05_0: 21.6517119661049
## log_srel_apap_Rawlins1977_1_0: 16.0027176603338
## log_srel_apap_Rawlins1977_2_0: 12.0075810603921
## log_srel_apap_sul_Chan1997_14_0: 142.487667782381
## log_srel_apap_sul_Chiew2010_56_0: 11.7517960003797
## log_srel_apap_sul_Critchley2005_14_0: 55.5656002676397
Some sigmas of the data are NA, try to recover an estimate for the sigmas from the other sigmas.
Take out some outliers for the fitting. Those are: study, name, time, reason 1. Chan1997, apap_cys, 5400, don’t know what went wrong with this one, but it just doesn’t fit reasonably in the time course 3. Chan1997, apap, 18000, The sigma is a few orders of magnitude lower. Mirjam said it might be that in this point less people were measured (eg 2) and they had nearly the same value. Then of course, sigma would be very smal
# data <- data %>%
# filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
# filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
# {.}
data %>% select(-n) %>% filter(study %>% str_detect("Chiew"))%>% as.datalist() %>% plotData()
data_with_errors <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name"), blather = T)
data_with_errors %>% select(study, name, D_apap, s0, srel) %>% unique()
## study name D_apap s0 srel
## 1 Chan1997 apap 1.4 -14.395145 -4.516431
## 13 Chan1997 apap_sul 1.4 -15.377176 -2.609451
## 26 Chan1997 apap_glu 1.4 -14.970003 -1.709818
## 39 Chan1997 apap_cys 1.4 -19.186945 -13.879472
## 51 Chiew2010 apap 5.6 -24.503629 -1.232431
## 65 Chiew2010 apap_glu 5.6 -9.710964 -2.303673
## 79 Chiew2010 apap_sul 5.6 -10.576758 -4.220755
## 93 Critchley2005 apap 1.4 -11.298969 -16.644666
## 108 Critchley2005 apap_sul 1.4 -16.328137 -2.802784
## 123 Critchley2005 apap_glu 1.4 -15.826025 -2.457868
## 138 Critchley2005 apap_cys 1.4 -18.851664 -18.687321
## 153 Rawlins1977 apap 0 -13.826728 -2.924789
## 163 Rawlins1977 apap 1 -28.577775 -1.968079
## 171 Rawlins1977 apap 2 -27.097492 -1.002046
## 179 Rawlins1977 apap 0.5 -27.374147 -2.115781
What’s the problem with the Chiew-dataset? The error model doesn’t work out, somehow
data_with_errors %>%
filter(study %>% str_detect("Chiew"), name %in% "apap") %>%
ggplot(aes(x=value)) +
geom_point(aes(y=sigmaLS^2*(n), color = log(time))) +
geom_line(aes(y=sigma^2*n)) +
geom_ribbon(aes(ymin=cbLower95, ymax=cbUpper95), alpha=.3) +
geom_ribbon(aes(ymin=cbLower68, ymax=cbUpper68), alpha=.3) +
ylab("variance") +
facet_wrap(~condidnt, scales = "free") +
scale_y_log10() +
theme_dMod() +
scale_color_continuous( low = "#98f5ff", high = "#4c4cdb")
myplot <- data_with_errors %>%
rename(sigma_fitted = sigma) %>%
left_join(data) %>%
gather("which_sig", "sigma", sigma, sigma_fitted) %>%
mutate(condition = paste0(study, D_apap)) %>%
ggplot(aes(x = time,y = log10(sigma))) +
geom_line(aes(color = condition, linetype = which_sig)) +
# geom_point(aes(color = condition, shape = which_sig))+
facet_wrap("name", scales = "free")
## Joining, by = c("study", "n", "time", "name", "value", "D_apap", "Ave_apap")
myplot %>% plotly::ggplotly()
# myplot
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters_e1 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters_e1 <- pars[!(names(pars)%in%c(free_parameters_e1,names(f)[1]))] %>% names
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_e1 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_e1 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_e1 <- c(free_parameters_e1, scale_parameters_e1)
# error_model_e1 <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_e1 <- setdiff(getSymbols(error_model_e1), names(error_model_e1)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_e1] <- paste0("exp(log", free_parameters_e1, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_e1, x = scale_parameters_e1, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_e1)
scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_e1, x = error_parameters_e1, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_e1)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_e1 <- NULL
for(i in 1:length(p_list)) { p_e1 <<- p_e1 + p_list[[i]]}
g_e1 <- Y(observables_e1, x)#, parameters = c(free_parameters_e1, scale_parameters_e1))
# err_e1 <- Y(error_model_e1, g_e1)
obj_e1 <- normL2(mydatalist, (g_e1*x*p_e1))#, errmodel = err_e1)
pouter_e1 <- rep(0, length(getParameters(obj_e1))) %>% set_names(getParameters(obj_e1))
# pouter_e1[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
#job_e1 <- runbg({myfit <- mstrust(objfun = obj_e1, center = pouter_e1, sd = 3, studyname = "methacetin", cores = 24, fits = 200); myfit}, machine = "ruprecht1", filename = "job_e1", input = global_env_without(c("mypred", "job", "myfit")))
# save(job_e1, file = "job_e1.rda")
# job_e1$check()
# load("job_e1_1_result.RData")
# myfit_e1 <- .runbgOutput
# load("job_e1.rda")
# myfit_e1 <- job_e1$get()
# save(myfit_e1, file = "myfit_e1.rda")
# job_e1$purge()
load("myfit_e1.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit_e1 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_e1 %>% as.parframe())#+scale_y_log10()
mypred_e1 <- (g_e1*x*p_e1)(mytimes*4, myfit_e1 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot_e1 <- plotCombined(mypred_e1, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_e1)
# myplot_e1
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 135 1151.996 TRUE 69 -18.255493 -2.247950
## 2 70 1151.996 FALSE 100 -25.864555 -2.246878
## 3 102 1151.996 TRUE 87 -20.902541 -2.247412
## 4 171 1151.996 TRUE 61 -15.814734 -2.245022
## 5 191 1152.002 TRUE 60 -7.374898 -2.255883
## 6 81 1152.021 FALSE 100 -5.978911 -2.280757
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -4.439263 -9.981200 -7.521667 -7.664488
## 2 -4.439394 -9.980696 -7.521636 -7.664391
## 3 -4.439279 -9.981359 -7.521708 -7.664394
## 4 -4.439563 -9.979954 -7.521437 -7.664217
## 5 -4.482247 -9.973992 -7.522330 -7.663926
## 6 -4.629608 -9.955660 -7.524306 -7.662562
## log_scale_apap_Chan1997_14_0 log_scale_apap_glu_Chan1997_14_0
## 1 0.6599173 12.8364130
## 2 0.6598512 20.4463668
## 3 0.6599152 15.4839438
## 4 0.6597140 10.3980699
## 5 0.6598690 1.9483154
## 6 0.6599017 0.5291009
## log_scale_apap_sul_Chan1997_14_0 log_scale_apap_cys_Chan1997_14_0
## 1 -1.0771287 -0.4991120
## 2 -1.0770584 -0.4996198
## 3 -1.0770897 -0.4989502
## 4 -1.0770215 -0.5003547
## 5 -1.0342097 -0.5064348
## 6 -0.8868909 -0.5250863
## log_scale_apap_Chiew2010_56_0 log_scale_apap_glu_Chiew2010_56_0
## 1 0.4747969 14.360289
## 2 0.4746995 21.970031
## 3 0.4747278 17.007724
## 4 0.4745103 11.921363
## 5 0.4764708 3.475552
## 6 0.4813080 2.066233
## log_scale_apap_sul_Chiew2010_56_0 log_scale_apap_Critchley2005_14_0
## 1 -0.8302279 0.4550537
## 2 -0.8301517 0.4549816
## 3 -0.8301895 0.4550346
## 4 -0.8301000 0.4548393
## 5 -0.7857834 0.4549590
## 6 -0.6342932 0.4548708
## log_scale_apap_glu_Critchley2005_14_0
## 1 13.2497621
## 2 20.8597371
## 3 15.8973054
## 4 10.8114742
## 5 2.3620052
## 6 0.9437402
## log_scale_apap_sul_Critchley2005_14_0
## 1 -1.373527
## 2 -1.373459
## 3 -1.373493
## 4 -1.373425
## 5 -1.330606
## 6 -1.183280
## log_scale_apap_cys_Critchley2005_14_0 log_scale_apap_Rawlins1977_0_1
## 1 -0.7385315 0.1580270
## 2 -0.7390379 0.1579934
## 3 -0.7383715 0.1580312
## 4 -0.7397729 0.1579200
## 5 -0.7458288 0.1579958
## 6 -0.7644113 0.1579937
## log_scale_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 0.3682622 0.1562957
## 2 0.3682065 0.1562248
## 3 0.3682818 0.1562731
## 4 0.3680777 0.1560860
## 5 0.3680259 0.1564571
## 6 0.3675088 0.1570966
## log_scale_apap_Rawlins1977_05_0
## 1 0.1238156
## 2 0.1237872
## 3 0.1238273
## 4 0.1236985
## 5 0.1235803
## 6 0.1229841
## logAPAPCYS_HLM_CL: 4.6261539108133e-05
## logAPAPCYS_Km: 0.000541229613400164
## logAPAPGLU_HLM_CL: 1.17961442872092e-08
## logAPAPGLU_Km: 0.105615516447963
## logAPAPSUL_HLM_CL: 0.0118046378599577
## logKa_apap: 0.000469196966461811
## log_scale_apap_Chan1997_14_0: 1.93463224874753
## log_scale_apap_Chiew2010_56_0: 1.60768769704826
## log_scale_apap_Critchley2005_14_0: 1.57625801184958
## log_scale_apap_cys_Chan1997_14_0: 0.607069478303199
## log_scale_apap_cys_Critchley2005_14_0: 0.47781508039879
## log_scale_apap_glu_Chan1997_14_0: 375649.964718959
## log_scale_apap_glu_Chiew2010_56_0: 1724226.50632285
## log_scale_apap_glu_Critchley2005_14_0: 567934.927599969
## log_scale_apap_Rawlins1977_0_1: 1.1711977984144
## log_scale_apap_Rawlins1977_05_0: 1.13180719602812
## log_scale_apap_Rawlins1977_1_0: 1.44522093867267
## log_scale_apap_Rawlins1977_2_0: 1.16917192088704
## log_scale_apap_sul_Chan1997_14_0: 0.340572019805132
## log_scale_apap_sul_Chiew2010_56_0: 0.435949924368047
## log_scale_apap_sul_Critchley2005_14_0: 0.253212324969791
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters_e2 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters_e2 <- pars[!(names(pars)%in%c(free_parameters_e2,names(f)[1]))] %>% names
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_e2 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_e2 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_e2 <- c(free_parameters_e2, scale_parameters_e2)
# error_model_e2 <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_e2 <- setdiff(getSymbols(error_model_e2), names(error_model_e2)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_e2] <- paste0("exp(log", free_parameters_e2, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_e2, x = scale_parameters_e2, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_e2)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_e2, x = error_parameters_e2, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_e2)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_e2 <- NULL
for(i in 1:length(p_list)) { p_e2 <<- p_e2 + p_list[[i]]}
g_e2 <- Y(observables_e2, x)#, parameters = c(free_parameters_e2, scale_parameters_e2))
# err_e2 <- Y(error_model_e2, g_e2)
obj_e2 <- normL2(mydatalist, (g_e2*x*p_e2))#, errmodel = err_e2)
pouter_e2 <- rep(0, length(getParameters(obj_e2))) %>% set_names(getParameters(obj_e2))
# pouter_e2[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# job_e2 <- runbg({myfit <- mstrust(objfun = obj_e2, center = pouter_e2, sd = 5, studyname = "methacetin", cores = 20, fits = 400); myfit}, machine = "ruprecht2", filename = "job_e2", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_e2, file = "job_e2.rda")
# job_e2$check()
# load("job_e2_1_result.RData")
# myfit_e2 <- .runbgOutput
# load("job_e1.rda")
# myfit_e1 <- job_e1$get()
# save(myfit_e2, file = "myfit_e2.rda")
# job_e1$purge()
load("myfit_e2.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit_e2 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_e2 %>% as.parframe())#+scale_y_log10()
mypred_e2 <- (g_e2*x*p_e2)(mytimes*4, myfit_e2 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot_e2 <- plotCombined(mypred_e2, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_e2)
# myplot_e2
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 379 2402.825 TRUE 56 -6.195743 -3.728993
## 2 125 2402.825 TRUE 46 -6.195869 -3.729162
## 3 147 2402.826 TRUE 68 -6.195203 -3.728064
## 4 224 2402.826 TRUE 54 -6.195879 -3.729148
## 5 374 2402.826 TRUE 39 -6.195582 -3.728692
## 6 28 2402.826 TRUE 38 -6.195853 -3.729130
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -5.922382 -10.65531 -7.895205 -7.156523
## 2 -5.922387 -10.65532 -7.895313 -7.156831
## 3 -5.922448 -10.65529 -7.894827 -7.156478
## 4 -5.922413 -10.65534 -7.895534 -7.156525
## 5 -5.922399 -10.65530 -7.895033 -7.156939
## 6 -5.922391 -10.65534 -7.895523 -7.156850
## logAPAPCYS_HLM_CL: 2.3575301587201e-05
## logAPAPCYS_Km: 0.000372525594735151
## logAPAPGLU_HLM_CL: 0.00203808917141075
## logAPAPGLU_Km: 0.0240169974041635
## logAPAPSUL_HLM_CL: 0.00267881249032963
## logKa_apap: 0.000779761067836627
The dynamics of e.g. apap_glu are too slow. Test two hypotheses:
This model also comprises some rate constants to make some reactions go faster
# Multiply rates from "ve" to "lu" with additional free parameters
reactions_to_lung <- reactions
reactions_to_lung[reactions_to_lung$from %>% str_detect("(1*Ave_apap)"),3] %<>% paste0(" * (k_to_lung_", c("sul", "", "cys", "glu"),")")
reactions_to_lung
el_to_lung <- eqnlist()
for(i in 1:nrow(reactions_to_lung)) el_to_lung <- addReaction(el_to_lung, reactions_to_lung$from[i], reactions_to_lung$to[i], reactions_to_lung$rate[i], reactions_to_lung$description[i])
# Convert to "eqnvec", which is basically a named vector of the ODEs and the names denote the states
f_to_lung <- el_to_lung %>% as.eqnvec()
parameters <- all_pars[getParameters(x)] %>% names()
free_parameters_to_lung <- parameters[parameters %>% sapply(. %>% str_detect(c("CL", "Km", "Ka_apap$")) %>% any) ]
# myodemodel_to_lung <- odemodel(f_to_lung, modelname = "methacetin_to_lung", fixed = setdiff(parameters, c(free_parameters_to_lung, "Ave_apap"))) # Include "Ave_apap" because of a bug: You need the sensitivities to at least one initial value
save(myodemodel_to_lung, file = "methacetin_to_lung.rda")
load("methacetin_to_lung.rda")
x_to_lung <- Xs(myodemodel_to_lung) # make prediction function
loadDLL(x_to_lung)
# get the only the parameters needed for x
pars_to_lung <- rep(1,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))
pars <- c(all_pars, pars_to_lung)[getParameters(x_to_lung)]
free_parameters_CLrenal <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap",
"CLrenal_apap", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_apap_sul"
)
fixed_parameters_CLrenal <- pars[!(names(pars)%in%c(free_parameters_CLrenal,names(f)[1]))] %>% names
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_CLrenal <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_CLrenal <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_CLrenal <- c(free_parameters_CLrenal, scale_parameters_CLrenal)
# error_model_CLrenal <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_CLrenal <- setdiff(getSymbols(error_model_CLrenal), names(error_model_CLrenal)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_CLrenal] <- paste0("exp(log", free_parameters_CLrenal, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_CLrenal, x = scale_parameters_CLrenal, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_CLrenal)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_CLrenal, x = error_parameters_CLrenal, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_CLrenal)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_CLrenal <- NULL
for(i in 1:length(p_list)) { p_CLrenal <<- p_CLrenal + p_list[[i]]}
g_CLrenal <- Y(observables_CLrenal, x_to_lung)#, parameters = c(free_parameters_CLrenal, scale_parameters_CLrenal))
# err_CLrenal <- Y(error_model_CLrenal, g_CLrenal)
obj_CLrenal <- normL2(mydatalist, (g_CLrenal*x_to_lung*p_CLrenal))#, errmodel = err_CLrenal)
pouter_CLrenal <- rep(0, length(getParameters(obj_CLrenal))) %>% set_names(getParameters(obj_CLrenal))
pouter_CLrenal[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
pouter_CLrenal[names(pouter_CLrenal) %>% str_detect("CLr")] <- -4
pouter_CLrenal[1:2] <- c(1.5,5)
# obj_CLrenal(pouter_CLrenal)
# job_CLrenal <- runbg({myfit <- mstrust(objfun = obj_CLrenal, center = pouter_CLrenal, sd = 3, studyname = "methacetin", cores = 16, fits = 100, iterlim = 150); myfit}, machine = "ruprecht2", filename = "job_CLrenal", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_CLrenal, file = "job_CLrenal.rda")
# job_CLrenal$check()
# myfit_CLrenal <- job_CLrenal$get()$ruprecht2
# # save(myfit_CLrenal, file = "myfit_CLrenal.rda")
# # job_CLrenal$purge()
load("myfit_CLrenal.rda")
loadDLL(x_to_lung)
## The following local files were dynamically loaded: methacetin_to_lung.so, methacetin_to_lung_s.so
myfit_CLrenal %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_CLrenal %>% as.parframe())+scale_y_log10()
mypouter <- myfit_CLrenal %>% as.parframe() %>% as.parvec
# mypouter[1:2] <- mypouter[1:2]+1
mypred_CLrenal <- (g_CLrenal*x_to_lung*p_CLrenal)(mytimes*4, mypouter, deriv = F)
# myplot_CLrenal <- plotCombined(mypred_CLrenal, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_CLrenal <- plotCombined(mypred_CLrenal, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_CLrenal)
# myplot_CLrenal
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 22 1679.527 TRUE 55 -6.532773 -4.069886
## 2 7 1679.527 TRUE 46 -6.533129 -4.070574
## 3 60 1679.527 TRUE 59 -6.532886 -4.069996
## 4 31 1679.527 TRUE 60 -6.532866 -4.069961
## 5 70 1679.527 TRUE 55 -6.532908 -4.070019
## 6 71 1679.527 TRUE 66 -6.532916 -4.070109
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logCLrenal_apap_sul
## 1 -5.629237 -10.64046 -6.378762 -5.099622
## 2 -5.629241 -10.64043 -6.374378 -5.099647
## 3 -5.629227 -10.64037 -6.382149 -5.099616
## 4 -5.629244 -10.64037 -6.383058 -5.099634
## 5 -5.629232 -10.64037 -6.383474 -5.099626
## 6 -5.629243 -10.64041 -6.380751 -5.099644
## logKa_apap logCLrenal_apap logCLrenal_apap_cys logCLrenal_apap_glu
## 1 -7.220531 -32.01661 -5.271552 -5.448181
## 2 -7.220536 -28.60979 -5.273667 -5.448045
## 3 -7.220632 -26.47175 -5.269753 -5.448253
## 4 -7.220556 -33.60964 -5.269310 -5.448251
## 5 -7.220625 -19.92413 -5.269103 -5.448267
## 6 -7.220526 -27.18079 -5.270518 -5.448182
## logAPAPCYS_HLM_CL: 2.39279181355532e-05
## logAPAPCYS_Km: 0.0016972225659718
## logAPAPGLU_HLM_CL: 0.00145496514549494
## logAPAPGLU_Km: 0.0170793409249032
## logAPAPSUL_HLM_CL: 0.00359131615388076
## logCLrenal_apap: 1.24556095423351e-14
## logCLrenal_apap_cys: 0.00513563209480641
## logCLrenal_apap_glu: 0.0043041255170805
## logCLrenal_apap_sul: 0.00609905356785048
## logKa_apap: 0.000731413761365896
Chiew2010 apap_sul looks good for the first time, even without scaling
load("methacetin_to_lung.rda")
x_to_lung <- Xs(myodemodel_to_lung) # make prediction function
loadDLL(x_to_lung)
# get the only the parameters needed for x
pars_to_lung <- rep(0,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))
pars <- c(all_pars, pars_to_lung)[getParameters(x)]
free_parameters_to_lung <- c("APAPGLU_HLM_CL", "APAPGLU_Km", "APAPSUL_HLM_CL", "APAPSUL_Km", "APAPCYS_HLM_CL", "APAPCYS_Km",
"Ka_apap",
paste0("k_to_lung_", c("sul", "", "cys", "glu"))
# "CLrenal_apap", "CLrenal_metc13", "CLrenal_apap_sul", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_co2c13"
)
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_to_lung <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_to_lung <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_to_lung <- c(free_parameters_to_lung, scale_parameters_to_lung)
# error_model_to_lung <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_to_lung <- setdiff(getSymbols(error_model_to_lung), names(error_model_to_lung)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_to_lung] <- paste0("exp(log", free_parameters_to_lung, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung, x = scale_parameters_to_lung, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_to_lung, x = error_parameters_to_lung, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_to_lung)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_to_lung <- NULL
for(i in 1:length(p_list)) { p_to_lung <<- p_to_lung + p_list[[i]]}
g_to_lung <- Y(observables_to_lung, x_to_lung)#, parameters = c(free_parameters_to_lung, scale_parameters_to_lung))
# err_to_lung <- Y(error_model_to_lung, g_to_lung)
obj_to_lung <- normL2(mydatalist, (g_to_lung*x_to_lung*p_to_lung))#, errmodel = err_to_lung)
pouter_to_lung <- rep(0, length(getParameters(obj_to_lung))) %>% set_names(getParameters(obj_to_lung))
pouter_to_lung[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
obj_to_lung(pouter_to_lung)
# job_to_lung <- runbg({myfit <- mstrust(objfun = obj_to_lung, center = pouter_to_lung, sd = 2, studyname = "methacetin", cores = 24, fits = 100, iterlim = 150); myfit}, machine = "ruprecht1", filename = "job_to_lung", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_to_lung, file = "job_to_lung.rda")
# job_to_lung$check()
# myfit_to_lung <- job_to_lung$get()$ruprecht1
# save(myfit_to_lung, file = "myfit_to_lung.rda")
# job_to_lung$purge()
load("myfit_to_lung.rda")
loadDLL(x_to_lung)
## The following local files were dynamically loaded: methacetin_to_lung.so, methacetin_to_lung_s.so
myfit_to_lung %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_to_lung %>% as.parframe())+scale_y_log10()
mypred_to_lung <- (g_to_lung*x_to_lung*p_to_lung)(mytimes*4, myfit_to_lung %>% as.parframe() %>% as.parvec(1) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
mypred_to_lung2 <- (g_to_lung*x_to_lung*p_to_lung)(mytimes*4, myfit_to_lung %>% as.parframe() %>% as.parvec(55) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F) %>% set_names(paste(names(.),"_2"))
# myplot_to_lung <- plotCombined(mypred_to_lung, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_to_lung <- plotCombined(c(mypred_to_lung,mypred_to_lung2), mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_to_lung)
# myplot_to_lung
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 70 1557.189 TRUE 142 1.882061 5.800306
## 2 34 1557.199 TRUE 148 1.296540 5.215847
## 3 41 1557.213 TRUE 142 1.652053 5.570782
## 4 4 1557.214 TRUE 150 1.600278 5.520080
## 5 54 1557.249 FALSE 150 1.315475 5.237448
## 6 95 1557.256 FALSE 150 1.215461 5.132848
## logAPAPSUL_HLM_CL logAPAPSUL_Km logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 1.3139749 5.074520 -7.348043 -6.117347
## 2 2.0624390 5.822052 -7.348164 -6.117514
## 3 1.0046156 4.764574 -7.348081 -6.117433
## 4 1.0356908 4.794728 -7.348121 -6.117488
## 5 0.7347650 4.491662 -7.348263 -6.117822
## 6 0.7286464 4.489584 -7.348001 -6.117180
## logk_to_lung_sul logKa_apap logk_to_lung_ logk_to_lung_cys
## 1 -0.1459383 -7.460770 -0.2039797 3.060832
## 2 -0.1448628 -7.460689 -0.2039490 3.060762
## 3 -0.1454166 -7.460765 -0.2039929 3.060800
## 4 -0.1444364 -7.460705 -0.2039785 3.060779
## 5 -0.1423281 -7.460575 -0.2039709 3.060692
## 6 -0.1464887 -7.460877 -0.2040240 3.060828
## logk_to_lung_glu
## 1 -0.8190345
## 2 -0.8201550
## 3 -0.8195732
## 4 -0.8207004
## 5 -0.8229671
## 6 -0.8182854
## logAPAPCYS_HLM_CL: 0.00064385123805241
## logAPAPCYS_Km: 0.0022042952730362
## logAPAPGLU_HLM_CL: 6.56702557689903
## logAPAPGLU_Km: 330.40080787627
## logAPAPSUL_HLM_CL: 3.72093453867929
## logAPAPSUL_Km: 159.895379325816
## logKa_apap: 0.000575213236929382
## logk_to_lung_: 0.815478948433962
## logk_to_lung_cys: 21.3453185074998
## logk_to_lung_glu: 0.44085707897448
## logk_to_lung_sul: 0.864211035946343
The plots look OK, but the values don’t make any sense.
x_to_lung <- Xs(myodemodel_to_lung) # make prediction function
loadDLL(x_to_lung)
# get the only the parameters needed for x
pars_to_lung_2 <- rep(0,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))
pars <- c(all_pars, pars_to_lung_2)[getParameters(x)]
free_parameters_to_lung_2 <- c("APAPGLU_HLM_CL", "APAPGLU_Km", "APAPSUL_HLM_CL", "APAPSUL_Km", "APAPCYS_HLM_CL", "APAPCYS_Km",
"Ka_apap",
paste0("k_to_lung_", c("sul", "", "cys", "glu"))
# "CLrenal_apap", "CLrenal_metc13", "CLrenal_apap_sul", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_co2c13"
)
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_to_lung_2 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_to_lung_2 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_to_lung_2 <- c(free_parameters_to_lung_2, scale_parameters_to_lung_2)
# error_model_to_lung_2 <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_to_lung_2 <- setdiff(getSymbols(error_model_to_lung_2), names(error_model_to_lung_2)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_to_lung_2] <- paste0("exp(log", free_parameters_to_lung_2, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung_2, x = scale_parameters_to_lung_2, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung_2)
scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
# scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_to_lung_2, x = error_parameters_to_lung_2, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_to_lung_2)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_to_lung_2 <- NULL
for(i in 1:length(p_list)) { p_to_lung_2 <<- p_to_lung_2 + p_list[[i]]}
g_to_lung_2 <- Y(observables_to_lung_2, x)#, parameters = c(free_parameters_to_lung_2, scale_parameters_to_lung_2))
# err_to_lung_2 <- Y(error_model_to_lung_2, g_to_lung_2)
obj_to_lung_2 <- normL2(mydatalist, (g_to_lung_2*x_to_lung*p_to_lung_2))#, errmodel = err_to_lung_2)
pouter_to_lung_2 <- rep(0, length(getParameters(obj_to_lung_2))) %>% set_names(getParameters(obj_to_lung_2))
pouter_to_lung_2[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# obj_to_lung_2(pouter_to_lung_2)
# job_to_lung_2 <- runbg({myfit <- mstrust(objfun = obj_to_lung_2, center = pouter_to_lung_2, sd = 2, studyname = "methacetin", cores = 24, fits = 100, iterlim = 150); myfit}, machine = "ruprecht1", filename = "job_to_lung_2", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_to_lung_2, file = "job_to_lung_2.rda")
# job_to_lung_2$check()
# myfit_to_lung_2 <- job_to_lung_2$get()$ruprecht1
# save(myfit_to_lung_2, file = "myfit_to_lung_2.rda")
# job_to_lung_2$purge()
load("myfit_to_lung_2.rda")
loadDLL(x_to_lung)
## The following local files were dynamically loaded: methacetin_to_lung.so, methacetin_to_lung_s.so
myfit_to_lung_2 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_to_lung_2 %>% as.parframe())+scale_y_log10()
mypred_to_lung_2 <- (g_to_lung_2*x_to_lung*p_to_lung_2)(mytimes*4, myfit_to_lung_2 %>% as.parframe() %>% as.parvec(35) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
# myplot_to_lung_2 <- plotCombined(mypred_to_lung_2, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_to_lung_2 <- plotCombined(mypred_to_lung_2, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_to_lung_2)
# myplot_to_lung_2
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 61 811.8280 FALSE 150 -6.998626 -17.50336
## 2 7 812.1984 FALSE 150 -7.019726 -16.97033
## 3 79 812.3622 FALSE 150 -7.031897 -16.94359
## 4 46 812.3956 FALSE 150 -7.033562 -16.92471
## 5 36 822.9495 FALSE 150 -7.920227 -17.47499
## 6 54 831.7163 FALSE 150 -9.112826 -16.08981
## logAPAPSUL_HLM_CL logAPAPSUL_Km logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 -4.647647 -7.067046 -8.005360 -25.35451
## 2 -4.677744 -7.095040 -8.017149 -24.29624
## 3 -4.699630 -7.149380 -8.015090 -24.25044
## 4 -4.703237 -7.161027 -8.014722 -24.18622
## 5 -4.908147 -8.078681 -7.579660 -23.58753
## 6 -5.018467 -8.363523 -7.356687 -23.10982
## logk_to_lung_sul logKa_apap logk_to_lung_ logk_to_lung_cys
## 1 -2.288579 -8.479607 -4.093847 -2.952564
## 2 -2.275876 -8.451360 -4.094141 -2.997499
## 3 -2.269219 -8.442444 -4.094653 -3.017023
## 4 -2.268317 -8.441194 -4.094752 -3.020282
## 5 -2.031370 -8.290519 -4.068860 -3.637812
## 6 -1.887041 -8.282598 -4.051058 -3.751590
## logk_to_lung_glu log_scale_apap_Chan1997_14_0
## 1 6.776023 2.035524
## 2 6.772015 1.963864
## 3 7.088785 1.967955
## 4 6.708936 1.971821
## 5 6.577030 2.122688
## 6 5.219719 2.144451
## log_scale_apap_glu_Chan1997_14_0 log_scale_apap_sul_Chan1997_14_0
## 1 6.495936 -2.067167
## 2 6.514004 -2.111232
## 3 6.844040 -2.125823
## 4 6.465985 -2.127919
## 5 7.263784 -2.351949
## 6 7.126559 -2.320297
## log_scale_apap_cys_Chan1997_14_0 log_scale_apap_Chiew2010_56_0
## 1 -4.984626 -0.1457100
## 2 -4.992699 -0.2872532
## 3 -5.003292 -0.3473295
## 4 -5.005029 -0.3556862
## 5 -5.625267 -1.0307136
## 6 -5.839268 -1.2249324
## log_scale_apap_glu_Chiew2010_56_0 log_scale_apap_sul_Chiew2010_56_0
## 1 8.738710 -2.787880
## 2 8.752954 -2.781700
## 3 9.079320 -2.776042
## 4 8.700771 -2.775228
## 5 9.406054 -2.621662
## 6 9.203844 -2.490344
## log_scale_apap_Critchley2005_14_0 log_scale_apap_glu_Critchley2005_14_0
## 1 1.914790 6.775920
## 2 1.843805 6.792860
## 3 1.848134 7.121887
## 4 1.852019 6.743747
## 5 2.008014 7.506618
## 6 2.029117 7.351299
## log_scale_apap_sul_Critchley2005_14_0
## 1 -2.517989
## 2 -2.561130
## 3 -2.575354
## 4 -2.577414
## 5 -2.786350
## 6 -2.748679
## log_scale_apap_cys_Critchley2005_14_0 log_scale_apap_Rawlins1977_0_1
## 1 -5.196222 -2.667753
## 2 -5.205343 -2.668026
## 3 -5.216401 -2.668411
## 4 -5.218228 -2.668487
## 5 -5.860658 -2.650692
## 6 -6.081074 -2.638401
## log_scale_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 2.251798 1.204153
## 2 2.148307 1.137746
## 3 2.142583 1.140532
## 4 2.145176 1.143994
## 5 2.125534 1.249054
## 6 2.134132 1.174536
## log_scale_apap_Rawlins1977_05_0
## 1 6.104588
## 2 5.154957
## 3 4.939859
## 4 4.913994
## 5 3.014949
## 6 2.848154
## logAPAPCYS_HLM_CL: 0.000333669227703504
## logAPAPCYS_Km: 9.74264737617257e-12
## logAPAPGLU_HLM_CL: 0.00091313568705377
## logAPAPGLU_Km: 2.50258813575314e-08
## logAPAPSUL_HLM_CL: 0.00958412454265644
## logAPAPSUL_Km: 0.000852748222445177
## logKa_apap: 0.000207660238495109
## logk_to_lung_: 0.016674967733492
## logk_to_lung_cys: 0.0522056534562965
## logk_to_lung_glu: 876.576074026939
## logk_to_lung_sul: 0.101410478546347
## log_scale_apap_Chan1997_14_0: 7.65626070083065
## log_scale_apap_Chiew2010_56_0: 0.864408357893709
## log_scale_apap_Critchley2005_14_0: 6.78551630713505
## log_scale_apap_cys_Chan1997_14_0: 0.00684233640219463
## log_scale_apap_cys_Critchley2005_14_0: 0.00553744269024556
## log_scale_apap_glu_Chan1997_14_0: 662.443889437685
## log_scale_apap_glu_Chiew2010_56_0: 6239.84412853667
## log_scale_apap_glu_Critchley2005_14_0: 876.485712326581
## log_scale_apap_Rawlins1977_0_1: 0.0694080270944431
## log_scale_apap_Rawlins1977_05_0: 447.908262977161
## log_scale_apap_Rawlins1977_1_0: 9.50480772104851
## log_scale_apap_Rawlins1977_2_0: 3.33393329597214
## log_scale_apap_sul_Chan1997_14_0: 0.126543722411801
## log_scale_apap_sul_Chiew2010_56_0: 0.0615515731003974
## log_scale_apap_sul_Critchley2005_14_0: 0.0806215896552064
Here I try to estimate the dosing. Dunno if this makes much sense, but it might be worth a try.
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters_est_input <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap" ,
"D_apap",
"CLrenal_apap", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_apap_sul"
#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters_est_input <- pars[!(names(pars)%in%c(free_parameters_est_input,names(f)[1]))] %>% names
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_est_input <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_est_input <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_est_input <- c(free_parameters_est_input, scale_parameters_est_input)
# error_model_est_input <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_est_input <- setdiff(getSymbols(error_model_est_input), names(error_model_est_input)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_est_input] <- paste0("exp(log", free_parameters_est_input, ")")
if(cond[1]==0) {
trafo["Ave_apap"] <- paste0("exp(logAve_apap_", rownames(conditions)[i], ")") %>% str_replace_all("\\.", "")
trafo["D_apap"] <- "0"
} else {
trafo["D_apap"] <- paste0("exp(logD_apap_", rownames(conditions)[i], ")") %>% str_replace_all("\\.", "")
trafo["Ave_apap"] <- 0
}
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung_2, x = scale_parameters_to_lung_2, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung_2)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_est_input, x = error_parameters_est_input, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_est_input)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_est_input <- NULL
for(i in 1:length(p_list)) { p_est_input <<- p_est_input + p_list[[i]]}
g_est_input <- Y(observables_est_input, x)#, parameters = c(free_parameters_est_input, scale_parameters_est_input))
# err_est_input <- Y(error_model_est_input, g_est_input)
obj_est_input <- normL2(mydatalist, (g_est_input*x*p_est_input))#, errmodel = err_est_input)
pouter_est_input <- rep(0, length(getParameters(obj_est_input))) %>% set_names(getParameters(obj_est_input))
pouter_est_input[names(myfit_CLrenal %>% as.parframe() %>% as.parvec())] <- myfit_CLrenal %>% as.parframe() %>% as.parvec()
pouter_est_input <- pouter_est_input[order(names(pouter_est_input))]
pouter_est_input[c(6,11:16)] <- log(c(1,1.4,5.6,1.4,0.5,1,2))
# p_est_input(pouter_est_input)
# job_est_input <- runbg({myfit <- mstrust(objfun = obj_est_input, center = pouter_est_input, sd = 5, studyname = "methacetin", cores = 16, fits = 200); myfit}, machine = "ruprecht2", filename = "job_est_input", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_est_input, file = "job_est_input.rda")
# job_est_input$check()
# (g_est_input*x*p_est_input)(times = mytimes, pouter_est_input, deriv = F) %>% plotPrediction(name %in% names(observables))
# load("job_est_input.rda")
# myfit_est_input <- job_est_input$get()$ruprecht2
# save(myfit_est_input, file = "myfit_est_input.rda")
# job_est_input$purge()
load("myfit_est_input.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit_est_input %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_est_input %>% as.parframe())+scale_y_log10()
mypred_est_input <- (g_est_input*x*p_est_input)(mytimes*4, myfit_est_input %>% as.parframe() %>% as.parvec(35) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
# myplot_est_input <- plotCombined(mypred_est_input, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_est_input <- plotCombined(mypred_est_input, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_est_input)
# myplot_est_input
## index value converged iterations logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 58 1122.299 FALSE 100 -10.48928 -4.957525
## 2 156 1122.300 TRUE 90 -10.49027 -4.961741
## 3 91 1142.864 TRUE 54 -10.48757 -5.127342
## 4 192 1142.864 TRUE 77 -10.48709 -5.124875
## 5 130 1142.864 TRUE 43 -10.48786 -5.128740
## 6 144 1142.864 TRUE 75 -10.48767 -5.127869
## logAPAPGLU_HLM_CL logAPAPGLU_Km logAPAPSUL_HLM_CL
## 1 7.615842 10.65704 -5.418030
## 2 6.535269 9.57645 -5.418026
## 3 10.713608 13.62289 -5.292063
## 4 18.348451 21.25771 -5.292033
## 5 12.851771 15.76106 -5.292060
## 6 9.264502 12.17378 -5.292044
## logAve_apap_Rawlins1977_0_1 logCLrenal_apap logCLrenal_apap_cys
## 1 0.10147172 -6.716203 -5.779922
## 2 0.10146728 -6.716345 -5.778142
## 3 0.02988996 -34.131247 -5.648947
## 4 0.02989775 -32.298981 -5.650080
## 5 0.02989460 -31.947337 -5.648331
## 6 0.02989872 -28.053280 -5.648702
## logCLrenal_apap_glu logCLrenal_apap_sul logD_apap_Chan1997_14_0
## 1 -5.444322 -4.590098 0.8099026
## 2 -5.444325 -4.590114 0.8099012
## 3 -5.331795 -4.452475 0.7546452
## 4 -5.331763 -4.452424 0.7546769
## 5 -5.331793 -4.452482 0.7546558
## 6 -5.331771 -4.452459 0.7546687
## logD_apap_Chiew2010_56_0 logD_apap_Critchley2005_14_0
## 1 2.563208 0.6401644
## 2 2.563184 0.6401447
## 3 2.467657 0.5631522
## 4 2.467679 0.5631632
## 5 2.467656 0.5631472
## 6 2.467666 0.5631538
## logD_apap_Rawlins1977_05_0 logD_apap_Rawlins1977_1_0
## 1 -0.8721721 0.12574694
## 2 -0.8721339 0.12575422
## 3 -0.9409273 0.05866417
## 4 -0.9409427 0.05867302
## 5 -0.9409016 0.05867707
## 6 -0.9408988 0.05868636
## logD_apap_Rawlins1977_2_0 logKa_apap
## 1 0.8141646 -7.855334
## 2 0.8141561 -7.855325
## 3 0.7695203 -7.811703
## 4 0.7695401 -7.811767
## 5 0.7695220 -7.811715
## 6 0.7695283 -7.811742
## logAPAPCYS_HLM_CL: 2.78331452999271e-05
## logAPAPCYS_Km: 0.00703030850292055
## logAPAPGLU_HLM_CL: 2030.10329039433
## logAPAPGLU_Km: 42490.6189263482
## logAPAPSUL_HLM_CL: 0.00443587568913188
## logAve_apap_Rawlins1977_0_1: 1.1067986135081
## logCLrenal_apap: 0.00121112777686759
## logCLrenal_apap_cys: 0.0030889577568275
## logCLrenal_apap_glu: 0.00432076935688847
## logCLrenal_apap_sul: 0.0101518605119084
## logD_apap_Chan1997_14_0: 2.24768909406063
## logD_apap_Chiew2010_56_0: 12.9773760893221
## logD_apap_Critchley2005_14_0: 1.896792740859
## logD_apap_Rawlins1977_05_0: 0.418042515914904
## logD_apap_Rawlins1977_1_0: 1.13399516723761
## logD_apap_Rawlins1977_2_0: 2.25728919295756
## logKa_apap: 0.000387678645051781
# Multiply rates from "ve" to "lu" with additional free parameters
reactions_to_lung_3 <- reactions
reactions_to_lung_3[reactions_to_lung_3$from %>% str_detect("(1*Ave_apap)"),3] %<>% paste0(" * (k_to_lung_", c("sul", "", "cys", "glu"),")")
reactions_to_lung_3
el_to_lung_3 <- eqnlist()
for(i in 1:nrow(reactions_to_lung_3)) el_to_lung_3 <- addReaction(el_to_lung_3, reactions_to_lung_3$from[i], reactions_to_lung_3$to[i], reactions_to_lung_3$rate[i], reactions_to_lung_3$description[i])
# Convert to "eqnvec", which is basically a named vector of the ODEs and the names denote the states
f_to_lung_3 <- el_to_lung_3 %>% as.eqnvec()
pars_to_lung_3 <- rep(1,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))
parameters <- c(all_pars, pars_to_lung_3)[getParameters(x)] %>% names()
free_parameters_to_lung_3 <- parameters[parameters %>% sapply(. %>% str_detect(c("_CL", "Km", "Ka_apap$", "k_to")) %>% any) ]
# myodemodel_to_lung_3 <- odemodel(f_to_lung_3, modelname = "methacetin_to_lung_3", fixed = setdiff(parameters, c(free_parameters_to_lung_3, "Ave_apap"))) # Include "Ave_apap" because of a bug: You need the sensitivities to at least one initial value
# save(myodemodel_to_lung_3, file = "methacetin_to_lung_3.rda")
load("methacetin_to_lung_3.rda")
x_to_lung_3 <- Xs(myodemodel_to_lung_3) # make prediction function
loadDLL(x_to_lung_3)
# get the only the parameters needed for x
pars <- c(all_pars, pars_to_lung_3)[getParameters(x)]
free_parameters_to_lung_3 <- c("APAPGLU_HLM_CL", "APAPGLU_Km", "APAPSUL_HLM_CL", "APAPSUL_Km", "APAPCYS_HLM_CL", "APAPCYS_Km",
"Ka_apap",
paste0("k_to_lung_", c("sul", "", "cys", "glu"))
# "CLrenal_apap", "CLrenal_metc13", "CLrenal_apap_sul", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_co2c13"
)
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_to_lung_3 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_to_lung_3 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_to_lung_3 <- c(free_parameters_to_lung_3, scale_parameters_to_lung_3)
# error_model_to_lung_3 <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_to_lung_3 <- setdiff(getSymbols(error_model_to_lung_3), names(error_model_to_lung_3)) %>% set_names(.,.)
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_to_lung_3] <- paste0("exp(log", free_parameters_to_lung_3, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung_3, x = scale_parameters_to_lung_3, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung_3)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_to_lung_3, x = error_parameters_to_lung_3, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_to_lung_3)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_to_lung_3 <- NULL
for(i in 1:length(p_list)) { p_to_lung_3 <<- p_to_lung_3 + p_list[[i]]}
g_to_lung_3 <- Y(observables_to_lung_3, x_to_lung_3)#, parameters = c(free_parameters_to_lung_3, scale_parameters_to_lung_3))
# err_to_lung_3 <- Y(error_model_to_lung_3, g_to_lung_3)
obj_to_lung_3 <- normL2(mydatalist, (g_to_lung_3*x_to_lung_3*p_to_lung_3))#, errmodel = err_to_lung_3)
getParameters(obj_to_lung_3)
myfit_to_lung %>% as.parframe() %>% as.parvec()
pouter_to_lung_3 <- myfit_to_lung %>% as.parframe() %>% as.parvec()
# job_to_lung_3 <- runbg({myfit <- mstrust(objfun = obj_to_lung_3, center = pouter_to_lung_3, sd = 4, studyname = "methacetin", cores = 24, fits = 1000, iterlim = 150); myfit}, machine = "ruprecht1", filename = "job_to_lung_3", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_to_lung_3, file = "job_to_lung_3.rda")
# job_to_lung_3$check()
# load("job_to_lung_3.rda")
# myfit_to_lung_3 <- job_to_lung_3$get()$ruprecht1
# save(myfit_to_lung_3, file = "myfit_to_lung_3.rda")
# job_to_lung_3$purge()
load("myfit_to_lung_3.rda")
loadDLL(x_to_lung_3)
## The following local files were dynamically loaded: methacetin_to_lung_3.so, methacetin_to_lung_3_s.so
myfit_to_lung_3 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_to_lung_3 %>% as.parframe())+scale_y_log10()
mypred_to_lung_3 <- (g_to_lung_3*x_to_lung_3*p_to_lung_3)(mytimes*4, myfit_to_lung_3 %>% as.parframe() %>% as.parvec(35) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
# myplot_to_lung_3 <- plotCombined(mypred_to_lung_3, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_to_lung_3 <- plotCombined(mypred_to_lung_3, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_to_lung_3)
# myplot_to_lung_3
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 598 1502.666 FALSE 150 5.222239 9.111693
## 2 585 1502.666 FALSE 150 6.010086 9.899766
## 3 230 1502.691 FALSE 150 7.145232 11.034067
## 4 74 1502.694 FALSE 150 3.539114 7.429065
## 5 901 1502.698 FALSE 150 3.013285 6.905417
## 6 852 1502.701 FALSE 150 2.820840 6.710850
## logAPAPSUL_HLM_CL logAPAPSUL_Km logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 4.908114 9.176861 -7.698137 -6.722977
## 2 3.700609 7.968713 -7.698307 -6.723640
## 3 1.294309 5.563838 -7.698262 -6.723571
## 4 2.076930 6.344381 -7.698409 -6.723928
## 5 3.797133 8.063931 -7.696645 -6.715827
## 6 5.819994 10.087409 -7.698370 -6.723744
## logk_to_lung_sul logKa_apap logk_to_lung_ logk_to_lung_cys
## 1 -0.4319789 -7.093422 -0.05287355 2.874983
## 2 -0.4313477 -7.093479 -0.05289277 2.874915
## 3 -0.4328742 -7.093376 -0.05285579 2.874963
## 4 -0.4307210 -7.093524 -0.05290436 2.874860
## 5 -0.4297375 -7.093220 -0.05277499 2.875319
## 6 -0.4306401 -7.093529 -0.05290571 2.874868
## logk_to_lung_glu
## 1 -0.7626255
## 2 -0.7628978
## 3 -0.7619719
## 4 -0.7632058
## 5 -0.7651345
## 6 -0.7632765
## logAPAPCYS_HLM_CL: 0.000453671805412169
## logAPAPCYS_Km: 0.00120295170864251
## logAPAPGLU_HLM_CL: 185.348656922024
## logAPAPGLU_Km: 9060.6175670097
## logAPAPSUL_HLM_CL: 135.383862495601
## logAPAPSUL_Km: 9670.74882908226
## logKa_apap: 0.000830550499923049
## logk_to_lung_: 0.948499944225183
## logk_to_lung_cys: 17.7251183888039
## logk_to_lung_glu: 0.466440188442027
## logk_to_lung_sul: 0.649223078801294
# save(list= ls(), file = "workspace.rda")